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Preface 


The  interaction  between  the  flight  control  system, 
structural  dynamics,  and  aerodynamics  of  aircraft  has  become  a 
major  concern  to  aircraft  designers.  There  has  been  considerable 
effort  to  develop  an  accurate  method  of  modelling  such 
interactions.  This  effort  develops  a  simple  method  to  create  an 
accurate  model  of  this  interaction.  This  model  can  then  be  used 
for  stability  and  control  analysis  including  the  effects  of 
structural  dynamics. 

I  wish  to  thank  my  advisor.  Dr.  Robert  A.  Calico  for  his 
invaluable  guidance  and  assistance  during  this  research  effort.  I 
would  also  like  to  thank  my  thesis  committee  members,  Dr.  Peter 
J.  Torvik  and  HaJ .  Lanson  Hudson  for  their  helpful  comments  and 
thorough  editing  of  this  document.  In  addition,  I  would  like  to 
thank  my  sponsor,  Dr.  Thomas  E.  Holl  for  his  guidance  on  applying 
aeroelastic  concepts  and  for  the  use  of  his  mathematical  model  of 
the  YF-iT.  Hy  thanks  also  go  to  Mr.  Maxwell  Blair  who  helped  me 
use  ADAM,  and  lLt.  William  Blake  who  helped  run  Digital  DATCOM. 
Finally,  I  wish  to  thank  my  wife  Colleen  for  Per  understanding 
and  support  during  this  work. 
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Abstract 


The  purpose  of  this  effort  was  to  provide  a  method  of 
developing  a  linear  model  of  an  elastic  aircraft.  The  model 
provides  the  capability  to  analyze  the  coupling  between  the  rigid 
and  elastic  motion  of  the  aircraft.  The  method  developed  in 
this  effort  obtains  stability  derivatives  directly  from  unsteady 
aerodynamic  forces.  This  results  in  a  state-space  model  whose 
states  are  just  the  normal  aircraft  states  and  rates,  the 
structural  coordinates  and  rates,  and  the  control  surface 
positions  and  rates.  Using  a  representation  of  the  YF-17  wind 
tunnel  flutter  model,  it  was  demonstrated  that  the  methodology 
developed  predicted  the  required  dynamics  to  make  this  a  viable 
method  of  modelling  rigid-body  and  flutter  behavior  of  the  model. 
Flutter  control  laws  were  designed  for  motion  about  an 
equilibrium  condition  represented  by  a  velocity  20X  above  the 
flutter  velocity.  Both  classical  and  modern  techniques  yielded 
acceptable  control  laws.  The  control  laws  were  also  analyzed  at 
off  design  conditions  to  check  robustness. 


MODELLING  OF  RIGID-BODY  AND  ELASTIC  AIRCRAFT  DYNAMICS 


FOR  FLIGHT  CONTROL  DEVELOPMENT 

I .  Introduction 

With  the  advent  of  high  gain  control  systems  on  high 
performance,  structurally  efficient  aircraft,  the 
interaction  between  aerodynamics,  structural  dynamics  and 
the  flight  control  system,  has  become  a  major  concern  to 
aircraft  designers.  This  interaction  has  been  termed 
aeroservoelasticity  (ASE).  Virtually  all  major  U.S.  fighter 
aircraft  in  use  today  have  encountered  this  phenomenon.  The 
F-15  has  encountered  several  ASE  problems  that  have 
influenced  both  ground  and  flight  tests  (10:Vol  I,  8-22). 

The  F/A-18  also  had  many  ASE  encounters,  primarily  in  ground 
tests  ( i 8: Vo  1  II,  205-211  ).  The  YF-16  and  YF-17  also  had  ASE 
encounters  (7:482).  Even  the  experimental  X-29  had  ASE 
problems,  although  not  in  flight.  The  ASE  interactions  are 
not  limited  to  new  aircraft.  The  ASE  phenomenon  has  been 
encountered  on  the  F-4  (18:Vol  I,  3  7;  7:482  )  and  many 
aircraft  that  were  in  test  phases  (18:Vol  II,  226-231; 
21:10).  This  interaction,  if  severe  enough,  can  limit  the 


performance  capabilities  of  an  aircraft.  Avoidance  of  such 


problems  calls  for  the  development  of  an  accurate  aircraft 
modelling  technique,  which  can  model  the  interactions 
correct  1 y . 


Background 

Typically,  the  three  phenomenon  involved  in  ASE,  have 
been  analyzed  using  separate  models  for  each.  The 
aerodynamic  models  are  usually  based  on  the  assumption  that 
the  aircraft  is  a  rigid-body.  An  accurate  representation  of 
the  aerodynamic  shape  of  the  aircraft  is  required  to  obtain 
accurate  aerodynamic  parameters.  The  flight  control  models 
are  typically  linearized  rigid-body  models,  using  stability 
derivatives  developed  for  a  rigid  aircraft.  These  linear 
flight  control  models  are  then  represented  by  state-space 
models  or  through  transfer  functions.  The  usual  way  to 
account  for  flexibility  is  to  add  an  elastic  correction 
factor  to  the  stability  derivatives,  which  does  not  account 
for  the  dynamics  of  the  structure.  The  structural  dynamics 
model,  on  the  other  hand,  is  concerned  mainly  in  a  correct 
representation  of  the  elastic  structure  to  predict  such 
things  as  flutter  and  divergence.  Typically,  each  of  the 
three  models  leads  to  an  analysis  which  is  never 
communicated  to  the  other  two.  This  lack  of  a  unified  model 


resulted  in  the  problems  that  were  noted  previously. 

Correcting  the  unwanted  ASE  interactions  after  the 
aircraft  is  built  and  tested  is  the  most  common  procedure 
for  dealing  with  ASE  problems.  The  typical  "quick  fixes" 
usually  involve  the  control  system.  Such  fixes  have 
included,  adding  filters,  reducing  control  system  gains, 
relocating  sensors  (18;  2i:8),  introducing  flutter  placards, 
and  in  the  case  of  the  X-29,  limiting  the  flight  envelope 
until  a  better  control  system  could  be  designed.  The 
process  of  separate  modelling,  and  then  fixing  the  problems 
after  the  aircraft  is  built  limits  the  potential  of  modern 
aircraft  design.  A  unified  approach  to  aircraft  modelling 
can  be  accomplished,  with  all  three  of  the  phenomenon 
involved  in  ASE  being  accounted  for. 

The  history  of  developing  a  unified  model  for  aircraft 
design  goes  back  well  over  twenty  years.  Bisplinghoff  and 
Ashley  (1:450-486)  and  Milne  (13)  were  some  of  the  first  to 
develop  the  equations  of  motion  for  an  elastic  aircraft. 

The  implementation  of  these  equations  was  beyond  the 
capabilities  of  computers  of  that  day.  In  1974,  the  first 
computer  program  that  addressed  the  ASE  phenomenon  was 
FLEXSTAB  (5).  FLEXSTAB,  developed  by  Boeing,  in  cooperation 
with  NASA  and  the  Flight  Dynamics  Laboratory  (FDL),  combined 
aerodynamic,  elastic  structure  and  control  system  models 
into  one  computer  program  for  analysis  and  design  purposes 


(5).  FLEXSTAB  was  a  large  and  cumbersome  program,  and  as  new 
design  and  analysis  techniques  were  developed,  they  could 
not  be  easily  incorporated  into  FLEXSTAB.  The  desire  to 
develop  a  unified  model  was  continued  throughout  the  70’s. 
Etkin  (6),  Schwanz  (22;  23),  Taylor  and  Woodcock  (26), 

Warren  (29),  Rodden  (19),  NASA  engineers  (18),  and  many 
others  (4;  11;  25)  all  proposed  approaches  to  develop  ASE 

models.  However,  it  became  apparent  that  even  by  1984  there 
was  no  consensus  on  the  best  way  to  develop  a  unified  model 
(18)  . 

Recently,  the  Flight  Dynamics  Laboratory  has  developed 
an  in-house  computer  program  called  ADAM  (Analog  and  Digital 
Aeroservoe 1 ast ic  Method).  It  has  the  capability  to  combine 
unsteady  aerodynamics,  multi-input  multi-output  control 
systems  and  a  structural  dynamics  model  into  one  analysis 
package  (14:  1  )  . 

ADAM  does  have  some  limitations.  The  first  limitation 
is  that  typical  stability  derivatives  that  flight  controls 
engineers  use  are  not  directly  available  (14:8).  A  second 
difficulty  is  that  ADAM  does  not  include  control  surface 
states  and  thus  cannot  predict  control  surface  ins tab i 1  1 1 i es 
(2:  Vol  1,  12).  In  addition,  ADAM  has  numerical  difficulties 

inherent  in  the  complexity  of  the  modelling  process  (14:8). 

A  method  is  needed  to  solve  these  problems  in  ADAM  to  make 
it  a  true  ASE  analysis  tool. 


The  purpose  of  this  thesis  is  to  develop  a  method  for 


use  in  concert  with  ADAM,  in  the  analysis  of  ASE  problems. 
This  method  includes  prediction  of  stability  derivatives, 
addition  control  surface  dynamics,  and  reducing  model 
comp  1  ex ity . 

I 

I 

Scope 


This  effort  will  be  limited  to  the  development  of  a 
linear  time- invariant  state-space  model  of  an  elastic 
aircraft.  This  effort  will  examine  only  longitudinal 
motion,  although  the  theory  can  be  easily  extended  to 
lateral  directional  motion.  The  work  is  also  limited  to 


continuous  time  models  and  control  systems. 


This  thesis  will  develop  the  equations  of  motion  of  an 
elastic  aircraft,  using  Lagrange’s  equations.  The  resulting 
equations  will  then  be  linearized  about  an  equilibrium 
point,  and  the  associated  non-linear  Aerodynamics  will  also 
be  linearized.  The  Kinematic  coupling  between  the 
rigid-body  and  elastic  motion  will  be  eliminated.  Thus  the 
only  coupling  between  the  rigid-body  and  elastic  motion  will 
be  through  the  aerodynamics.  The  resulting  second  order 
equations  will  then  be  transformed  into  a  first  order 
state-space  model  for  the  elastic  aircraft.  A  computer 
program  will  then  be  developed  which  will  automate  this 
procedure.  The  methodology  will  be  demonstrated  with  a 
representation  of  the  YF-17  wind  tunnel  flutter  model. 


II.  Theoretical  Development 


It  is  desired  to  develop  the  equations  of  motion  of  an 
elastic  aircraft,  in  the  form  of  Eqs  i  and  2 


X  = [A] x+ [B] u 

(1  ) 

Y--  [C ]  x 

(2  ) 

where 

x  -  n  vector  of  states  of  the  aircraft 

u  -  m  vector  of  inputs  to  the  aircraft 

Y  ~  p  vector  of  outputs  of  the  aircraft 

[A]  -  n  by  n  plant  matrix 

[B]  -  n  by  m  input  matrix 

[C]  -  p  by  n  output  matrix 

There  are  three  general  ways  that  the  equations  of 
motion  for  an  elastic  aircraft  can  be  determined.  The  first 
way,  and  probably  the  most  common  way  is  to  use  classical 
Newtonian  dynamics.  EtKin  (6:122-145)  developed  the 
equations  of  motion  for  a  rigid  aircraft  using  Newtonian 
dynamics.  Milne  (13:4),  Taylor  (26:22),  and  Bisplinghoff 
and  Ashley  (1:450-466)  developed  the  equations  of  motion  for 
an  elastic  aircraft  using  this  method.  The  second  method  is 
to  use  Hamilton’s  equations  to  derive  the  equations  of 
motion  (3:1684-1687).  The  final  method  of  deriving  the 
equations  of  motion  uses  Lagrange’s  equations,  as  Schwanz 
had  done  (22).  This  is  the  method  used  in  this  development, 
since  it  is  the  most  direct  and  easy  to  use.  Lagrange’s 


equations  contain  all  the  information  needed  to  derive  the 
equations  of  motion. 


Lagrange » s  Equations 

Lagrange's  equations  are  based  on  energy  principles. 

If  the  aircraft's  total  potential  and  kinetic  energy  can  be 
described,  in  reference  to  an  inertial  frame,  then  the 
equations  of  motion  of  that  aircraft  can  be  derived.  Stated 
in  a  general  form,  Lagrange's  equations  for  a  holonomic 
system  are 


d/dt ( dL/da l - dL/dg= Q 

(3) 

L  =  T-V 

(4) 

where 

T  -  total  kinetic  energy, 

V  -  total  potential  energy, 

3  -  generalized  coordinates, 

Q  -  generalized  forces,  and 
dL/d^:  (dL/dfit  .  .  .  dL/d-nnJT 

Therefore  in  order  to  use  Lagrange's  equations,  the  kinetic 
and  potential  energy  of  the  elastic  aircraft  must  be  found. 

Consider  an  elastic  body,  as  shown  in  Figure  1.  The 
inertial  axis  origin,  Oj,  is  a  fixed  point  attached  to  the 
earth.  It  will  be  assumed  that  the  earth  is  flat  and 


non-rotating.  In  general  the  flat  non-rotating  earth 


assumption  does  not  have  to  be  made,  however,  it  is 
reasonable  to  do  so  for  aircraft  motion  in  subsonic  and  low 
supersonic  flight  (6:148).  This  assumption  will  allow  the 
notation  to  be  kept  simple,  and  simplifies  the  dynamics. 

The  position  vector  from  this  inertial  origin  to  any 
arbitrary  point  P  on  the  body, R^  is  given  by: 


SpzBo*£o*~  (5) 

where 

Rq  -  position  vector  of  the  body  frame  orig 
to  the  inertial  frame  origin, 
rQ  -  position  vector  of  the  undeformed 

position  of  P  to  the  body  axis  origin, 

6  -  elastic  deformation  position  vector  of 
point  P  from  the  undeformed  position. 


From  this  the  kinetic  and  potential  energies  can  be 
obtained . 


in 

and 


Kinetic  Energy .  The  kinetic  energy  of  the  body  is 
def ined  as 


Trl/a/pRp1 . RpJd  V 
V 


(6) 


where 

p  -  density 
V  -  volume, 

Rp1  -  velocity  of  point  P  with  respect  to  Oj 
as  seen  by  an  observer  in  the  inertial 
reference  frame 


can  be  represented  by 
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£pI=ioI+dI/dt<r0+d) 


(7  ) 


which  can  be  expanded 


RpI=ioI+AB+QB/Ix(r0+d)  (6) 

where 

qB/I  -  angular  velocity  of  the  body  axis  with 
respect  to  the  inertial  axis 
jB-  velocity  of  the  deformed  point  P  with  respect 
to  the  body  axis  frame 

Substituting  this  back  into  Eq  6  yields 


T=i/2jp{V0+<5B+QB/Ix(r0+dH  .  {Vo+dn+O8/^  (r0+d  ) }  d  V'  (9  ) 
V 


where,  RqI=Yo.  Expanding  Eq  9  term  by  term  gives 


T=l/2j,P¥o-Yodvr+1/2«fp¥o  •iBdK+i/2j'pV0  .  (Q^xfro  +  d  )  ]<1V 

V  v  V 

+  l/2j‘pdB.V0dV+l/2j'pdB.dBdV'+  1/2^00®.  (fiB/Ix(r0id)|d7 

V  v  V 

+  l/2j'p{QB/Ix(r0+d) )  .y0d  V+  l/2j'piOB/Ix  (r0+d  ) )  ,aBdV 

V  V 

♦  l/2j'p{0B/Ix  (  r0+d  )  J  .{QB/Ix(r0+d))dV'  (10) 

V 


The  total  mass  of  the  aircraft  is  defined  as 


(il ) 


Z/rJ'pdV' 

V 

Combining  term  in  Eq  10,  and  expanding  the  remaining 
products  yields 

T=1/2MV0  •Vo+J'pXo  V+/pV0  .  (0B/Ixr0  JdK+J'pVo  .  (QB/Ixa  )d  V 
V  V  V 

+  i/2j'p<$B.6Bd  V+XpgB/1  .  (r0xdB)d  V+J'pqB/1  .  (dxdBjdK 

V  V  V 

+  1/2  Sp  (QP/lxr0  )  .  (qB/^Pq  )d  w/p(QB/Ixr0  )  .  (QB/IX4  )d/ 

V  V 

+  J*P (qB/Ix^ ) . (qB/Ix^ )d y  (12) 

V 

Since  Vq  and  0B/I  d0  not  depend  on  the  position  within 
the  volume,  they  can  be  pulled  outside  the  integrals.  Now 
the  following  assumptions  are  made  to  help  simplify  Eq  12. 
First,  it  is  assumed  that  the  body  axis  origin  is  at  the 
center  of  mass  of  the  equilibrium  configuration  of  the 
aircraft.  This  implies  that 


0  =  vfprod  V  (13) 

V 

The  second  major  assumption  is  that  the  aircraft  body  axis 
system  is  the  mean  axis  system,  which  Milne  states  as  the 
axis  n.  .  .at  every  point,  the  linear  and  angular  momentum 
of  the  relative  motion  with  respect  to  the  body  axis  is 


wwm  mm  uj^xRy<s^ir 


identically  zero'*  (13:5).  This  results  in  the  following 


0  =  J'pdd  V-  J*pr  0xdd  vtfpdBd  Vz /proXdBd  V  (14) 

V  ~  V  ~  V  V 

This  reduces  the  coupling  between  the  overall  motion  and  the 
elastic  deformation  (13:27),  thus  reducing  the  Kinetic 
energy  equation  (Eq  12)  to 

T  =  1  / 2W0  .  v0+  1  /2  JpdB  .  6  Bd  y+  j*py0  .  (  qB/  I x 6  )  d  V 

V  V 


+  l/2j‘p(QB/Ixr0  )  .  (flB/^xro  )d^+j‘p(QB/Ixr0  )  .  (QB/IXd  )<1V 
V  V~ 

+  J‘p(QB/Ix<5)  .  (QB/Ixd)dVr  (15) 

V 

Noting  that  the  components  of  the  cross  product  of  two 
vectors  expressed  in  the  same  orthogonal  reference  frame, 
say  axb  =  c,  may  be  obtained  from 

[a] b=c  (16) 

where  [a]  is  a  skew  symmetric  matrix,  defined  by 


•j\ri 

I 


[a] : r  0  -a3 

I  a3  0 
L-a2  a! 


Eq  15  now  becomes 


“al  I 
OJ 
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I 

{ 


T=i/ZfN0.V0+i/zSp±B.±BdV*i/ZQB/1  .  [IR]  .O®/1 

V 

+  Sf>QB/1  .  (6x±B)4V+Qb/1  .  [IC]  .Q®/1 

V 

+1/2Q®/1 . [IE] .Q®/1  (17 ) 

where 

[I  ]  -  rigid  body  inertia  matrix 

[I*]  -  inertial  coupling  matrix  between  rigid  and 
v  elastic  motion 

[lg]  -  elastic  motion  inertia  matrix 

It  is  at  this  point  that  Knowledge  about  the 

equilibrium  conditions  is  necessary  to  simplify  Eq  17.  Since 

it  is  desired  to  develop  a  linear  model,  perturbation 

equations  about  an  equilibrium  condition  will  be  developed. 

The  equilibrium  condition  in  this  development  will  be 

straight  and  level  flight  of  an  undeformed  airframe,  thus 

the  equilibrium  body  rotation  rates  (p  ,q  ,r  )  are  zero. 

e  e  e 

Thus,  Q®/1  is  now  Just  a  vector  of  the  perturbation 
rotation  rates.  Both  the  coupling  and  elastic  inertial 
terms  of  Eq  17  are  functions  of  <5  ,  a  perturbation 
motion.  Thus  the  last  three  terms  of  Eq  17  are  third, 
third,  and  fourth  order  respectively  of  the  perturbation 
motion,  and  thus  will  be  ignored.  Thus,  Eq  17  becomes 

T=l/2JW0.V0+l/2j'pd®.d®dV=H/2g®/I  .  [IR]  .O®/1 

V 


(18) 


much  simpler  than  the  original  equation  (E q  10).  The  only 
term  that  involves  the  elastic  motion  is  the  middle  term  in 
Eq  18.  It  is  at  this  point  that  some  assumptions  about  the 
elastic  deformations  must  be  made.  In  general,  the  elastic 
deformations  of  a  body  can  be  described  by  a  set  of  n 
coupled  second  order  ordinary  differential  equations, 
representing  the  free  vibration  of  the  body.  These  coupled 
equations  can  be  decoupled  by  use  of  a  linear  transf ormat ion 
(12:143-144).  Therefore,  in  this  development,  it  will  be 
assumed  that  a  set  of  uncoupled  orthogonal  modes  are 
available.  This  assumption  will  allow  diagonal  mass  and 
stiffness  matrices  to  be  developed  for  the  equations  of 
motion,  thus  allowing  easier  decoupling  of  the  elastic  and 
rigid  body  equations  of  motion.  The  elastic  deformations 
are  defined  as 

n 

i  =  <X,  y,  2  )4i  (t )  ( 19  ) 

1  =  1 

where, 

Ti(x,  y,  z)  -  mode  shape  of  the  ith  mode,  and 
5l(t)  -  time  dependent  motion  of  the  mode 

Since  the  modes  are  orthogonal,  JVi?j=0  for  i^J 
(12:  143).  Now  defining  the  generalized  mass  as 


(20) 


M-i  =Xp«r  i2<i  v 

V 


and  substituting  Eqs  19  and  20  into  Eq  18  yields 


T=l/2iW0.V0+l/2Eui4i2+l/2QB/1  .  [IR]  . Q®/  1 
i  =  1 


(21  ) 


This  is  the  equation  that  is  required  for  the  Kinetic 


energy 


Potential  Energy .  The  potential  energy  of  an  elastic 


structure  can  be  represented  by 


V=-l/2ZUi^i2Wi2 
lr  l 


(22) 


where 


wj2  -  the  squared  natural  frequencies  of  the  various  modes 


The  same  assumption  about  orthogonality  that  were  made 


earlier  is  also  made  here  (12:188).  The  potential  energy  due 


to  the  Earth's  gravitation  will  be  dealt  with  as  an  external 


force  acting  on  the  aircraft. 


Eqs  21  and  22  can  now  be  substituted  into  Eq  4,  the 


equation  for  L,  and  Lagrange’s  equations  formed,  and  noting 


that  the  rigid-body  terms  do  not  appear,  gives 
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ft 


ft  V- 
•aV 

Pi< 

i 

ivr, 


k: 

i 


m 

1 


& 


(23  ) 


n 

EtUiii  +  HiWi2*!  )=Q 
1  =  1 


This  can  be  represented  in  matrix  notation  as 


where 

li}  -  ith  generalized  mass 
-  ith  modal  frequency 
Si  -  ith  generalized  modal  coordinate 
F|  -  ith  generalized  force 

In  this  formulation,  the  control  surface  motions  will 
be  treated  as  extra  degrees  of  freedom  in  the  elastic 
motion.  In  order  to  correctly  develop  the  remaining 
equations  for  rigid  body  motion,  the  definitions  of  linear 
and  angular  momentum  must  be  used.  The  linear  momentum  may 
be  found  from 

P=dT/dV0  (25) 

Then  the  equations  of  motion  for  rigid  translation  become 
FRrdVdttPJrdVdttdT/dVo) 

=dB/dt(dT/dV0)+OB/Ix(dT/dV0 )  (26) 

Using  the  typical  formulation,  the  aircraft  equations  of 


TU*|  r 


I  r  0  -pi  I V I  =  I  Y+Fy-Mgs in<pcos0 
L-q  p  oj  LwJ  LZ+Fz-MgCOS*COS0j 


(27  ) 


where 


M 

u 

V 

w 

p 

q 

r 

X 

Y 
Z 

F„ 


mass 

velocity  in 
velocity  in 
velocity  in 
roll  rate 
pitch  rate 
yaw  rate 
aerodynamic 
aerodynamic 
aerodynamic 
-  thrust  in 


the 

the 

the 


direction 

direction 

direction 


force 
force 
f  orec 
the  x 


in  the  x 
in  the  y 
in  the  z 
direction 


direction 

direction 

direction 


Fy  -  thrust  in  the  y  direction 


F i  -  thrust  in  the  z  direction 

of  gravity 


g  -  acceleration 
6  -  pitch  angle 
V  -  bank  angle 


with  the  following  equations  for  kinematics  and  trajectory 


Ji 


1  0  -sin0T 

0  cosy  sintcosol 

0 

.0  -sin*  costcosoj 

y. 

c0cy 

c0sy 

-sOl 

’x' 

s*s 0cy-c*sy  sfs0sy+ctcy 

S*C0 

y 

.c*s0cy-s*sy  cvsOsy-s+cy 

C*C0J 

.z. 

(28) 


(29) 


where 


c  -  cosine 
s  -  sine 
y  -  yaw  angle 


A  similar  method  can  be  used  to  develop  the  equations 
of  motion  for  the  attitude  motion.  The  angular  momentum  is 
given  by 


-  18  - 


HrdT/dO®/1 

The  equations  for  rigid  body  rotation  become 
M0=dI/dt (H)=dI/dt (dT/bO®/1 ) 

=d®/dt (dT/dG®/1 )+0®/Ix(dT/dQB/1 ) 


(31  ) 


These  equations  (Eqs  26  and  31)  are  sometimes  termed 
Lagrange's  equations  in  quasi-  coordinates.  Eq  31  in  matrix 
notation  becomes 


[*xx  -Ixy  _Ix2l 
"Jxy  Jyy  'Jyz  9 
-IXj  “lyi  IzzJ lrJ 

'  0  -r  ql  r  Ixx  -IXy  -Ixzl  TPl  rL] 
♦  r  0  -P  -Ixy  1 1 yy  -Iyz  U  M 

L-q  P  Oj  L-IX2  -Iyx  IxzJjrJ  LhJ 


(32) 


where 


I  XXi  I  x  y .  IX2i  I  yy  ^  yz*  1 2  2  mass  moment 
or  product  of  inertia 
L  -  aerodynamic  moment  about  the  x  axis 

M  -  aerodynamic  moment  about  the  y  axis 

N  -  aerodynamic  moment  about  the  z  axis 


Eqs  24,  27,  28,  29,  and  32  comprise  the  desired  equations  of 
motion.  The  forces  in  Eqs  24,  27,  and  32  are  now  function 

of  both  the  rigid  and  elastic  motion. 


Longitudinal  Equations  of  Motion 


The  desired  equations  of  motion,  Eqs  i  and  2,  are 


it 
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linear  in  nature.  The  above  equations  (Eqs  27,  28,  29,  and 

32)  still  contain  non-linear  terms  which  need  to  be  further 
linearized.  In  this  development,  however,  only  the 
longitudinal  equations  will  be  developed.  Therefore,  only 
those  equations  appropriate  to  the  longitudinal  motion  will 
be  addressed.  The  equilibrium  condition  is  again  straight 
and  level  flight,  but  now  using  longitudinal  motion, 
equations  reduce  to  the  following  scalar  equations 


«U+>#qw: -X+Fx-/#gs  ine 

(33  ) 

Mir - Mu  u  =  -  Z  +  F  2  -  Afg  c  o  s  0 

(34) 

I yyd 

(35) 

q  =  0 

(36) 

u=xcos6-zs  in0 

(37  ) 

w  =  xs  in8+  zcose 

(38) 

The  attitude  equations  are  uncoupled  from  the  trajectory 

equations  (Eqs  37  and  38),  and  thus  the  trajectory  equations 

will  be  disregarded  in  this  development.  Now  it  will  be 

assumed  that  each  variable  will  be  the  sum  of  its’ 

equilibrium  value  and  perturbed  value  (i.e.  u  =  u  tu1  ),  and 

e 

the  derivatives  of  the  equilibrium  values  are  zero.  It  will 
also  be  assumed  that  the  mass  and  inertia  remain  essentially 
constant.  If  the  thrust  is  to  be  used  as  an  input  it  can  be 
used  as  a  control  input.  Substituting  these  into  Eqs  33-38 


(39) 


Ml*+Af(q'  )(we+w»  )=-(Xe+X'  )  +  Fx-/fgs  in  (0e+0  *  ) 

q’  )(ue  +  u»  )  =  -(Ze  +  Z*  )  +  Fz-/fgcos  (0e+e»  )  (40) 

Iyyq'=Me+M*  (41) 

< q *  )  =  e  (42) 

u,=x,co*(0e*0' )-i » s in (8e+0 »  )  (43  ) 

u» =x» sin (0^+0* )  +  z 1  cos ( 0 e+ 0 '  )  (44) 


Now  taking  all  vhe  terms  that  strictly  deal  with  equilibrium 
conditions  yields  Eqs  45-47  while  the  remaining  terms  result 
in  the  perturbation  equations,  Eqs  48-53. 


0  =  -Xe+Fx-/fgs  in0e 

(45) 

O  =  -Ze*Fz-/fgcos0e 

(46) 

Me  =  0 

(47  ) 

/ni+JVlWeq'+w*  q*  )  =  -X '  -Mgs  in9 * 

(48) 

/Hs+l‘f(ueq,  +  u,q'  )  =  — Z  * 

(49) 

Iyyq*  =M» 

(50) 

q  ’  =  0  ’ 

(51  ) 

c 

II 

X- 

(52  ) 

W'  =2  ’ 

(53  ) 

The  equilibrium  equations  and  the  perturbation  equations  are 
identical  to  those  developed  in  Taylor  (26:39-44).  Assuming 
that  the  terms  where  products  of  perturbation  quantities 
exists  are  small  results  in  the  linear  set  of  equations  in 


K 


>y\ 

CSS: 
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Eq  54,  where  the  primed  superscript  is  dropped. 


'MOO 

fvl 

■ 

OHO 

* 

♦ 

,0  0  1  yy_ 

Lq. 

0  0  Mmt 
0  0  Mje 
0  0  0 


(54) 


The  equations  can  be  combined  with  the  elastic  equations  of 
motion  to  yield  the  following  matrix  equation: 


(55) 


These  equations  are  identical  to  those  of  Schwanz  (23:54). 
The  linear  aerodynamics  which  drive  these  equations  must  now 
be  determined  to  complete  the  model  of  the  flexible 
aircraft.  From  this  point  on,  only  the  short  period  motion 
will  be  considered.  Thus  the  first  equation  in  Eq  55  will 


not  be  considered  and  there  will  not  be  an  x  state. 


Linear  Aerodynamics 


Since  the  equations  of  motion  of  the  aircraft  have  been 
linearized,  the  aerodynamics  associated  with  this  model  must 
also  be  linearized.  Etkin  (6:159)  and  Taylor  (26:45)  have 
both  shown  that  the  aerodynamics  can  be  approximated 
linearly.  Ignoring  symmetric  and  antisymmetric  motion 
interactions,  and  neglecting  most  higher  order  aerodynamic 
derivatives,  the  longitudinal  aerodynamic  forces  are  of  the 
f o I  towing  form 


m 

HrMuU’+Mqq'+Maa+MQa+ElH^jdj+M^jdj+Mjjdj ) 

J  "  ^ 

n 

♦  S(M€  +  +  £i)  (56) 

1:1  1  1  1 

where, 

MyrdM/du,  etc. 

In  general  these  derivatives  are  not  directly  available 

from  the  unsteady  aerodynamic  methods  such  as  the  Doublet 

Lattice  (8),  which  is  used  in  flutter  prediction  work. 

Rodden  and  Giesing  (20)  have  shown  however,  that  stability 

derivatives  can  be  obtained  from  such  methods.  In 

developing  unsteady  aerodynamic  forces  and  moments,  it  is 

assumed  that  the  motion  is  purely  oscillatory  and  of  the 
iwt 

form  e  The  resulting  generalized  forces  therefore 


contain  real  and  imaginary  parts.  The  aerodynamics  are 


determined  at  a  constant  altitude  and  Mach  number  for 
various  reduced  frequencies  (K:wb/V,  where  b  is  a 
reference  semi-chord  and  V  is  the  velocity).  This  results 
in  a  set  of  forces  for  each  degree  of  freedom  for  each 
reduced  frequency. 

One  way  to  approximate  the  frequency  dependent  is  to 
use  a  Pade  fit  as  a  function  of  reduced  frequency  (14:3-4). 
In  general  the  equations  of  motion  can  be  represented  by  the 
second  order  equations 


[M)a+eC]s*[K]3=-(pV2s/2)[Q<lC)J3  (57) 

where 

[M]  -  generalized  mass  matrix 
[C]  -  generalized  damping  matrix 
[K]  -  generalized  stiffness  matrix 
[Q ( K ) ]  -  generalized  aerodynamic  matrix 
3  -  generalized  coordinate  vector 
V  -  velocity 
S  -  reference  area 


The  [Q(K)]  matrix  is  a  matrix  of  complex  coefficients  based 
on  simple  harmonic  motion.  In  order  to  solve  the  flutter 
problem  in  a  classical  way,  many  different  values  of  k  and 
the  associated  [Q(k)J  matrices  are  required.  An  alternate 
technique  is  to  curve  fit  each  coefficient  in  the  [Q ( K ) ) 
matrix  with  respect  to  K,  using  Pade  polynomials.  The 


polynomials  are  of  the  form 


(58) 


.»  *•!  I  i  »1  |  tj‘>  I  i. 


La*i  m'k  i’t  t 


.t ,  JL.ii.fr, 


Qij(k)  =  tH0+Hi  (iK)  +  N2(ik)2  +  N3  (ilC)3  +  N4(ik)4l/ 

(1  +  Nq+Dj  (ik)+D2(ik)2) 


for  a  fourth  order  over  second  order  fit.  Then  substituting 


bs/V=ik  into  Eq  58  yields 


Qij (s )={Mo+Ni (bs/v)+M2<b*/v>2+N3 (bs/V)3+N4<bs/V)4)/ 

{1+Nq+Dj (bs/V)+D2(bs/V)2) 


(59) 


Taking  the  Laplace  transform  of  the  right  side  of  Eq  57,  and 
equating  terms  of  like  power  in  s,  results  in  the  following 
equation 


1  [I]sil+[a3)s3+[a2]s2+ta1)s+[a0n3(s)=0  (60) 


As  can  be  seen,  this  method,  increases  the  order  of  the 
system  by  the  multiple  associated  with  the  order  of  the 
denominator  of  the  Pade  fit.  The  system  order  is  therefore 
altered  unless  the  Pade  denominator  is  a  constant.  This 
makes  it  extremely  difficult  to  form  a  state-space  model  of 
the  aircraft  in  which  all  the  states  relate  directly  to 
physical  quantities.  Another  method,  suggested  by  Rodden 
and  Giesing  (20)  allows  the  state  variables  to  be  explicitly 
stated,  and  may  also  reduce  the  order  of  the  model.  This  is 
the  method  that  will  be  used  in  this  study. 
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In  order  to  derive  the  equations  necessary  to  find  the 


derivatives  some  assumptions  must  be  made.  In  the  analysis 


of  the  unsteady  aerodynamic  forces,  the  rigid-body  motions 


are  assumed  to  be  pitch  (q)  and  plunge  (h).  In  aircraft 


dynamics,  plunge  is  typically  replaced  with  angle  of  attack 


(a).  Plunge  is  pure  vertical  motion,  and  is  related 


directly  to  w  by  (h  =  w)  if  measured  from  the  center  of 


mass.  If  a  is  small,  which  can  be  assumed  since 


perturbation  motion  is  being  discussed,  then;  a=w/V=h/V. 


And  since  the  equilibrium  condition  assumes  level  flight 


then  a= e  for  the  equilibrium  conditions 


The  method  of  Rodden  and  Giesing  was  developed  to 


calculate  dynamic  stability  derivatives  from  unsteady 


aerodynamics  (20).  This  methodology  is  summarized  here, 


using  the  moment  equation  as  an  example.  First,  it  is 


assumed  that  stability  derivatives  can  be  represented  by  a 


Maclaurian  series 


M=M0  +  Maa+M(ja+Mga+Mqq+ 


(61  ) 


Oscillatory  theory  cannot  predict  the  M0  term,  and 


therefore  it  will  be  omitted.  Next  it  is  assumed  that  the 


pitching  and  plunging  motion  are  oscillatory  and  are  of  the 


form  a:0:eoeiut  and  h  =  h0elll,t,  and  thus  for 


pitching  motion 


SL- 


M=Real (Me iwt ) 


(62  ) 


Hence  the  complex  pitching  amplitude  M  due  to  pitching 
motion  can  be  represented  by  the  series  expansion  to  first 
order 


H=e0tM<x+i«(M<j,+Mq)] 

Similarly  for  plunging  motion  with  a  =  li/V 


(63a) 


H=h0/2c  ( iwMaVMfl,)  (63b) 

The  terms  of  Eqs  63  can  now  be  directly  related  to  the 
appropriate  terms  of  the  matrix  of  complex  frequency 
dependent  aerodynamic  coefficients  from  the  Doublet  Lattice 
method.  According  to  Rodden  and  Giesing,  the  steady  dynamic 
stability  derivatives  are  defined  by  the  limiting  values  as 
k  approaches  0,  so  in  general,  a  small  value  of  k  should  be 
used.  For  greater  detail  on  the  procedure  used  by  Rodden 
and  Giesing,  the  reader  is  referred  to  their  paper  (20). 

The  application  of  the  procedure  outlined  above  in  this 
effort  is  as  follows.  The  unsteady  aerodynamic  forces  were 
obtained  from  FASTOP  (30)  which  uses  Doublet  Lattice 
aerodynamics.  The  generalized  forces  for  a  reduced 
frequency  are 


-  27  - 


vV 


(64) 


QtJ=  (1/S2  )J\f<APj/q)  (hi/S  )dS 
S 

where 

Q  -  the  force  on  the  ith  mode  due  to  the  jth 
J  modal  deflection 
APj  -  pressure  on  the  jth  mode 
h  -  magnitude  of  the  ith  mode 
s 1 -  reference  length 
S  -  surface  area 
q  -  dynamic  pressure 

In  order  to  obtain  stability  derivatives  the  Q.  terms 

i  1 

must  be  normalized  by  both  the  ith  and  jth  mode  shape 
magnitudes  as  was  demonstrated  by  Noll  with 

lateral-directional  derivatives  (15).  To  develop  dimensional 
stability  derivatives  the  Cr^  terms  must  be  multiplied  by 
the  dynamic  pressure  (q),  and  the  reference  area  (S).  Also, 
any  Q.^  terms  which  are  related  to  moments  must  also  be 
multiplied  by  a  reference  length  (c  -  reference  chord). 

Thus  the  terms  in  Eqs  65-73  are  now  dimensional  quantities. 
Eqs  65  are  the  rigid'body  results  of  Rodden  and  Giestng  in 
dimensional  form 


Pitc 

h  and 

Plunge  Influence 

on  "Rigid  Forces" 

QZZ  = 

(  iu>Za 

-u>2Zd)  (h0/V)  (h0  )e 

iut 

( 65a  ) 

Qmz2 

( iu>Ma 

-u^Ma)  (h0/V)  (0o)e 

iuit 

(65b  ) 

qZM= 

1  i<*>  (  Zq+  Zq  )  1  ( 60  )  ( h0  ) 

eia>t 

( 65c  ) 

Qhm1 

IMa*  lUi(M6+Mq)  )  (  0  0  )  (0O  ) 

e  iut 

( 65d  ) 
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Sodden  and  Giesing's  method  for  rigid-body  motion  can  easily 


be  extended  further  to  obtain  the  elastic  and  control 
surface  influences  on  the  rigid  body  forces  and  visa  versa, 
along  with  the  interactions  among  themselves.  These  are 
shown  in  Eqs  66-73. 


Elastic  Inf luences  on  "Rigid  Forces" 

QZCj3  <z61+lwZ^1‘w2zl1 )(5i0)(ho>«iWt 

i  =  (MC  / ♦  1 1  )  ( e  o  )  e  11 wt 

Pitch  and  Plunge  Influence  on  "Elastic  Forces" 

^  iZ=  <  i«F^  ta-w2Fc  i<3, )  (h0/V )  ( 6  1q  )e  i“»t 

i«(Fgid*Fgiq))  (•0)Ui#)«iwt 


Elastic  Inf luence  on  "El ast ic  Forces" 

Qe  e  =<f£  £  +iwFe  £  -u^F £  s  )  ( 6 1  )  ( 6 j  )« 

i  J  i  J  i  “i  i 


iwt 


(66a) 

(66b) 

(67a) 

(67b) 

(68) 


Contro 1  Surface  Inf  1 uence  on  "Rigid  Forces" 

<*Z<J  =(Z<5.  +  iwz<J.-w2z3.  )(<5i  )(h0)eiwt  (69a) 

i  i  i  i  o 

<>Hdi=(M6.  +  iwM8i-w2M3i  )  <8^  )  (0O  )e  iwt  (69b) 

Pitch  and  Plunge  Inf luence  on  "Contro 1  Surface  Forces" 
Q<,iZ=(iwFdia-w2F<5i(5,)(h0/V)(alo)eiwt  (70a) 

G^M'fF^a*  iw(F<ji<3,+  F<ji<I)|  (0o  )  (<J )e  ita)t 


(70b) 


Contro 1  Surface  Inf  1 uence  on  "El ast 1 c  Forces " 

Qf  <5  =(F^  4  +  iwF^  a  -u>2Fc  3  )(3i  )Uj  )elwt  (71) 

sj  l  J  l  J  i  J  i  l  J 

Elastic  Inf luence  on  "Contro  1  Surf  ace  Forces" 

Q<5  *  +lwFd  £  -u^F^  £  )Ui  )(<5j  )e^  (72) 

J1  J1  J1  J  ^  ^  J 

Contro 1  Surf  ace  Inf luence  on  HContro 1  Surf  ace  Forces" 

Q<s  a  MFa  i  +iu>Fa  a  -uzFd  3  M<*i  )(«Sj  )elwt  (73) 

J  i  j  1  j  i  j  1  i  J 

In  FASTOP,  to  find  the  steady  derivatives,  R  is  set  to  zero, 
thus  only  the  real  terms  of  the  above  equations  appear.  The 
remaining  derivatives  are  found  at  a  small  value  of  R.  The 
expansion  was  Rept  to  second  order  which  results  in  no 
increase  in  system  order.  If  higher  order  dynamics  are 
needed,  they  can  be  easily  added  by  adding  higher  order 
stability  derivatives.  This  will,  of  course  increase  the 
order  of  the  system. 


Resu  1 1  ini 


luat  ions 


The  stability  derivatives  from  Eqs  65-73  can  now  be 
moved  to  the  left  hand  side  of  Eq  55,  thus  filling  the 
matrices  with  terms .  The  only  terms  left  on  the  right  hand 
side  of  Eq  55  are  the  generalized  forces  due  to  the  torques 
applied  to  the  control  surface  hinges  from  the  control 


inputs,  as  shown  in  Eq  74. 


It  is  shown  Kane  (9:81-83)  that  the  generalized  forces,  Q, 
due  to  the  control  torques  can  be  represented  by 


m 

C^rEdO/dTii  .Tj  (75) 

j  =  l 


where 


T 


j 


This  can 


torques  applied  to  the 

be  represented  by 


control 


surface  hinges 


Q=  [dQ/diu]  [Ti  •  •  •  Tm]T  (76) 

The  torque  can  be  related  to  the  input  signal  to  the  surface 
by 


TJ=n,jaJ 


where 


(77  ) 


m  -  the  generalized  mass  of  the  ith  control  surface 
6l-  is  the  command  signal  to  the  ith  control  surface 


Now  letting  8=  [8j  .  .  .<$m]  Eq  75  becomes 


Q= [E]u 


(78  ) 


where 


[EJrtdO/dai] [diag  mj] 


.■•icyt.ti.tt  j1  * 


Thus  Eq  74  becomes 


[X] x+ [Y]x+ [Z]x  = [E]u 


(79  ) 


If  the  matrix  [X]  can  be  inverted,  then  Eq  79  can  be 
rearranged  into  the  typical  state-space  representat ion 


g]=[-x-?z  -x-i.]g]*[-x-?Ja 


(80) 


Eq  80  is  the  desired  equation  in  the  form  of  Eq  1.  Eq  80  and 
the  procedure  for  obtaining  the  dimensional  stability 
derivatives  from  the  generalized  forces  (Eqs  65-73)  were 
programmed  into  a  computer  procedure  called  MAC  (Methodology 
for  Aeroe l ast ic lty  and  Controls).  Appendix  B  is  the  User’s 
Manual,  and  Appendix  C  are  the  subroutines  that  comprise 


32  - 


& 

K 


1 


III.  YF-17  Mode  1  Deve I opment 

In  order  to  demonstrate  the  methodology  developed  in 
the  previous  section,  a  simple  model  that  exhibited  flutter 
characterist ics  was  needed.  A  mathematical  representation 
of  the  YF-17  wind  tunnel  flutter  model  was  used.  This  model 
had  the  required  structural  properties  needed  to  validate 
the  methodology.  Second,  only  the  first  two  elastic  modes 
were  necessary  to  adequately  predict  flutter.  This  results 
in  a  simple  representation  of  the  model  dynamics. 

The  model  used  is  shown  in  Figure  2.  It  is  represented 
by  six  modes;  the  rigid  body  pitch  and  plunge  modes,  two 
structural  modes  (first  wing  bending  and  first  wing  torsion 
modes),  and  two  control  surface  modes  (the  trailing  edge 
flap  and  an  all  moveable  horizontal  tail).  This  model  will 
allow  full  testing  of  the  influences  desired.  Table  I 
contains  the  generalized  masses  and  structural  frequencies 
for  the  model.  The  model  developed  includes  an  elastic 
wing,  connected  by  a  rigid  extension  to  a  rigid  tail.  The 
mode  shapes  for  the  elastic  modes  were  obtained  from  Noll's 
work  (IS).  The  geometric  and  structural  data  were  entered 
into  FASTOP,  a  flutter  prediction  program  (30).  FASTOP  was 
used  to  determine  the  generalized  forces  for  reduced 


frequencies  of  0.0  and  0.02.  The  flutter  prediction 
capability  of  FASTOP,  which  uses  the  P-K  method  of  flutter 
prediction,  was  also  used  to  check  the  validity  of  this 
prediction  method.  The  flutter  prediction  capability  of 
ADAM  (2:Vol  1,  24-25),  which  uses  a  ’velocity  root  locus’ 

technique  to  predict  flutter,  was  also  used  as  a  check. 

Two  sensors  were  chosen  near  the  trailing  edge  flap,  at 
40.154  inches  along  the  span,  as  shown  in  Figure  2.  Each 
sensor  measures  trans 1  at ional  motion  at  that  wing  location. 
The  two  outputs  are  then  differenced  and  divided  by  the 
difference  between  the  two  sensors,  to  obtain  the  twist 
angle  at  that  span  station.  The  resulting  output  matrix, 

[C]  is  then  a  combination  of  the  two  structural  mode  shapes 
at  that  span  station.  Appendix  A  contains  the  output  matrix 
used.  The  sensor  and  actuator  dynamics  are  not  modelled  in 


this  development. 


>,»’»*»>»l»>l^  iU’i*.  »1.  .f. 
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IV.  Results  and  Discussion 


The  resulting  equations  (Eq  BO)  are  linearized 


approximat ions  to  the  actual  aircraft  dynamics.  Where  the 


aerodynamics  are  linear  with  respect  to  the  reduced 


frequencies,  this  model  should  successfully  capture  the 


essential  dynamics  and  thus  predict  the  low  frequency 


flutter,  with  the  added  benefit  of  predicting  rigid  body 


stability  derivatives.  The  resulting  state  space  model  will 


also  be  of  lower  order  than  that  obtained  using  aerodynamics 


modelled  using  the  Pade  approxlmant  method,  since  it  will 


not  contain  the  often  unnecessary  higher  order  aerodynamics 


associated  with  curve-fitting  the  aerodynamics  in  the 


frequency  domain. 


In  order  to  prove  that  this  method  is  a  useful  way  to 


model  elastic  aircraft  dynamics  for  flight  control  use, 


three  tasks  were  analyzed.  First,  the  ability  of  this 


method  to  predict  stability  derivatives,  and  thus,  the 


ability  to  predict  rigid  body  dynamics,  was  investigated. 


Next,  the  methodology's  ability  to  predict  flutter  and 


associated  structural  instabilities  was  investigated 


Finally,  the  model  was  used  to  develop  a  workable  flutter 


control  law. 
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STABILITY  DERIVATIVE  COMPARISON  FOR  THE  YF-17  MODEL 


Derivative 


Da  t com 


cLa 

4.305 

4.142 

cMa 

-0.6501 

-0.5038 

cLa 

2 . 623 

2.276 

cMa 

-4.699 

-7 .312 

CLq 

6.03 

8.145 

cMq 

-6.799 

-7 .474 

cLde 

1 .054 

1  .  153 

CM<5 

e 

-1 .886 

-1 .909 

CL,( 

0 . 1432 

0.233 

CMdf 

-0.0751 

-0.4366 

Stability  Derivative  Prediction 


The  capability  of  MAC  to  predict  rigid  body  stability 
derivatives  was  validated  using  the  YF-17  model.  The 
stability  derivatives  from  MAC  were  compared  to  the 
analytical  methods  of  Digital  Datcom  (28).  As  can  be  seen  in 
Table  II,  MAC  does  an  acceptable  Job  in  predicting  most 
stability  derivatives,  with  the  exception  of  the  flap 
derivatives.  If  there  are  better  sources  for  stability 
derivatives,  MAC  has  the  capability  to  incorporate  them  into 
the  anal ys is . 
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Flutter  Prediction 


The  capability  of  MAC  to  predict  flutter  was  also 
evaluated  with  the  YF-17  model.  The  results  from  MAC  were 
compared  to  two  different  flutter  prediction  methods.  The 
first  is  FASTOP,  which  uses  the  P-K  method  of  predicting 
flutter  speeds  and  frequencies  (30).  The  second  is  ADAM, 
which  uses  a  'velocity  root  locus'  technique  to  predict 
flutter  speeds  and  frequencies  (2:Vol  1,  24-25).  Figure  3  is 
a  plot  of  ADAM'S  root  locus.  MAC  uses  a  similar  technique 
using  Eq  80.  Figure  4  shows  the  plot  of  the  velocity  root 
locus  of  MAC  for  various  speeds.  As  can  be  seen  by  Figures 
3  and  4,  both  methods  are  in  very  close  agreement  in 
predicting  the  dynamics  of  the  structural  modes.  All  three 
methods  agree  very  well  in  predicting  the  flutter  speed  and 
frequency  (Table  III).  Using  FASTOP  as  a  basis,  it  can  be 
seen  that  both  ADAM  and  MAC  slightly  underpredict  the  speed 
of  flutter  while  slightly  overpredicting  the  frequency.  All 
the  predictions  are  within  5*  of  one  another. 

In  comparing  the  eigenvalues  of  MAC  with  ADAM  (Table 
IV)  it  is  seen  that  the  higher  order  aerodynamic  roots 
associated  with  the  Pade  fit  method  of  ADAM  do  not  seem  to 
be  Important  to  the  dynamics  involved.  This  is  similar  to 
the  result  that  Pasquini  reported  (17:31-32). 
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TABLE  III 


COMPARISON  OF  FLUTTER  PREDICTION  METHODS  FOR  THE  YF-17  MODEL 


Flutter 

Prediction 

Method 


Flutter 
Speed 
(ft/sec ) 


Flutter 
Frequency 
(Hertz ) 


FASTOP 

ADAM 

MAC 


6 . 10 
6.25 
6.18 


TABLE  IV 


EIGENVALUE  COMPARISON  OF  MAC  WITH  ADAM 


FOR  THE  458  FT/SEC  CASE 


Eigenvalue 


ADAM 


Rigid  Body 


Elastic 


Control  Surface 


0.000 

0.000 

-2 . 412+5. 586i 
-2.412-5.5861 
-5.890+34.331 
-5. 890-34. 33i 
3 .371+36.671 
3 . 371-36 . 67 1 
-0.518+314.41 
-0.518-314.41 
-12 .93+373 .3  1 
-12 .93-373 .3 i 


838+3 .3141 
838-3 .3141 
215+6. 084i 
215-6.0841 
901+34.261 
901-34. 26i 
485+36 . 66 i 
485-36 . 661 


Aerodynamic 


93 .92+ 184.7  i 
93 .92-184.7 i 
94.01+184.01 
94.01-184.01 
94.34+184.21 
94.34-184.21 
94.23+184.31 
94.23-184.21 
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Imaginary  Component 


Control  Design 

From  the  previous  two  sections  it  can  be  concluded  that 
Eq  80  accurately  captures  the  essential  dynamics  of  the 
YF-17  model.  The  model  from  MAC  was  then  used  to  develop  a 
flutter  control  system,  using  both  classical  and  modern 
control  techniques.  The  benefit  of  using  this  method  of 
developing  the  aircraft  equations  of  motion  is  that  the 
state-space  model  (Eq  80)  is  formed  from  the  second  order 
form  of  the  equations  of  motion  (Eq  79).  The  state  vector 
then,  is  comprised  of  the  generalized  coordinates  and  their 
rates.  In  contrast,  if  the  Pade  approximation  method  is 
used,  the  state-space  model  must  be  formed  from  the  transfer 
function  equation  (Eq  60).  Typically,  this  introduces 
additional  states,  and  yields  states  which  are  not  the 
generalized  coordinates  and  their  rates. 

Both  control  designs  were  based  on  the  YF-17  model  at 
1.2  times  the  flutter  speed  (458  ft/sec).  All  the 
eigenvalues  of  this  condition  are  in  Table  V. 

Classical  Design.  Transfer  functions  were  determined 
from  the  state-space  model  for  the  above  condition  using 
MATRIXx  (10).  The  transfer  function  from  the  sensor  to  the 


flap  (y/<Sf  )  was  then  used  to  develop  a  feedback 


TABLE  V 


EIGENVALUES  OF  THE  UNAUGMENTED  YF-17  AT  458  FT/SEC 


Eigenvalue  Frequency  (Hx ) 


0.000  0.000 

0.000 

“2.412+5. 5861  0.889 

-2.412-5.5861 

-5.890+34.331  5.464 

-5.890-34.331 

3. 371+36. 67i  5.836 

3.371-36.671 

-0.518+314.41  50.04 

-0.518-314.41 

-12.93+373.31  59.42 

-12.93-373.31 


law.  Using  root  locus  techniques  in  TOTAL  (27),  the  gain 
was  varied  until  the  unstable  structural  mode  became  stable 
while  keeping  the  remaining  structural  mode  and  the  short 
period  roots  stable.  The  gain  chosen  was  -  50.  Figure  5  is 
the  root  locus  for  the  y/de  transfer  function,  with 
twist  position  feedback.  As  can  be  seen,  the  unstable 
structural  mode  is  made  stable,  however,  the  very  lightly 
damped  control  surface  mode  goes  slightly  unstable.  This  is 
not  of  much  concern  since  the  control  surface  actually  has 
damping  due  to  frictional  forces,  which  are  not  accounted 
for  in  this  model.  The  resulting  roots  of  this  condition 
with  the  control  system  are  in  Table  VI. 


This  simple  control  system  was  then  tested  at  off 
design  conditions.  As  can  be  seen  from  the  ’velocity  root 


x -poles 


o- zeros 

A- root  locus  values  for  a  gain  of  -50 


Figure  5.  Root  Locus  of  y/&f  of  the  YF-17  Model 
at  458  ft/sec 


TABLE  VI 


EIGENVALUES  OF  THE  YF-17  MODEL  AT  458  FT/SEC 
USING  TWIST  POSITION  FEEDBACK 


Eigenvalue 


0.000 
0.000 

-2 .826+1 .5791 
-2.826-1 .5791 
-3  .003+35.541 
-3.003-35.541 
-0 .082+37 .441 
-0.082-37.441 
0.460+312.11 
-0.460-312.11 
-12.92+373.31 
-12.92-373.31 


Frequency  (Hz ) 


0.000 


0.251 


5.  S56 


5.959 


49.67 


59.42 


locus’  In  Figure  6,  the  control  system  keeps  the  structural 
modes  and  the  short  period  mode  stable  for  all  conditions 
above  250  ft/sec  up  to  the  design  condition,  thus  yielding 
an  acceptable  control  system,  if  it  is  implemented  after  the 
model  reaches  250  ft/sec. 


Modern  Design.  Again  using  the  1.2  flutter  speed 


I'M 


-6.  -4.  —  2. 


Real  Component  (1/sec) 


Figure  6. 


Velocity  Root  Locus  for  the  YF-17 
Model  with  Twist  Position  Feedback 


TABLE  VII 


EIGENVALUES  OF  THE  YF-1T  MODEL  AT  458  FT/SEC 


USING  8TATE- SPACE  FEEDBACK 


Eigenvalue 


Frequency  (Hi ) 


0.000 

8.41845.5941 
8 .418-5.5941 
5.889434.331 
5.889*34.331 
3.3T5438.6T1 
3.3T5-36.8T1 
10.004314.41 
10 .00*314.41 
18.9343T3.31 
18.93-3T3.3i 


0.890 


5.464 


5.836 


50.04 


59.48 


-  45  - 


MM 


HP 


The  feedback  Matrix  waa  then  used  at  the  off  deaifn 
condition! .  Thla  atate-apace  feedback  ayateM  kept  the 
atructural  and  abort  period  Modea  atable  for  all  the 
condition!  above  100  ft/aec,  up  to  the  deaifn  condition!, 
and  even  beyond  (Figure  T). 

Quite  often  the  atatea  are  not  readily  Meaaurable. 
However,  the  Mode  11  inf  technique  aa  outlined  in  thia  effort 
can  be  uaed  in  the  deaifn  of  an  eatiMator  or  observer.  Thia 
Model  can  be  uaefully  eMployed  in  a  Modern  control  deaifn. 


-6.  -4  -Z  0  Z.  4.  6 


Real  Component  (1/sec) 

Figure  7.  Velocity  Root  Locus  for  the  YF-17 
Model  Using  State-space  Feedback 


V.  Conclusions  and  Recommendations 

The  equations  of  motion  of  an  elastic  aircraft,  have 
been  developed  in  this  thesis,  including  the  effects  of 
control  surface  dynamics.  As  was  shown,  the  model  developed 
is  useful  for  stability  and  control  analysis,  including  the 
effects  of  structural  dynamics.  The  model  can  be  used  as  a 
basis  for  either  classical  or  modern  control  methods.  This 
model  also  has  merit  as  a  flutter  prediction  method,  at 
least  where  the  aerodynamics  are  well  behaved,  and  can  be 
linearly  approximated  at  low  reduced  frequencies. 

There  are  several  advantages  to  using  this  method. 
First,  stability  derivatives  can  be  predicted,  but  if  a 
better  source  of  stability  derivatives  is  Known,  they  can  be 
directly  incorporated  into  the  model.  A  second  advantage  is 
that  the  states  are  the  normal  aircraft  states  and  rates, 
plus  those  associated  with  the  structural  generalized 
coordinates  and  rates,  and  those  associated  with  the  control 
surface  rotations  and  rates.  This  makes  this  method 
directly  applicable  to  classical  and  modern  control  design 
procedures.  A  third  advantage  is  that  the  second  order 
approximation  used  in  this  thesis  results  in  a  lower  order 
model  than  the  methodology  used  in  ADAM.  Even  if  higher 
order  dynamics  are  needed,  they  can  be  easily  added  to  the 
model  as  pointed  out  by  Rodden  (20).  Finally,  this  method 


can  allow  the  prediction  of  control  surface  instabilities, 
and  unlike  ADAM  can  predict  rigid  body  notion. 

This  project  is  a  useful  design  tool  and  it  should 
continue  to  be  expanded.  The  use  of  other  dynanic  models 
should  be  used  to  further  validate  this  method,  including 
models  in  which  the  rigid  modes  couple  with  the  structural 
modes.  Second,  the  ability  to  extend  this  model  to  beyond 
second  order,  and  to  extend  it  to  lateral-directional 
motion.  Finally,  the  ability  to  account  for  sensor  and 
actuator  dynamics,  along  with  control  system  dynamics  needs 
to  be  addressed,  to  provide  a  better  dynamic  model. 


YF-17  State-Space  Model  at  1 . 2  V 


This  appendix  contains  the  A,  B,  and  C  matrices  for  the 
YF-17  model  at  20  percent  above  the  flutter  speed  (458 
ft/sec).  The  states  of  the  system  are  in  Table  Ai .  The  A 
matrix  is  contained  in  Table  A2 .  The  B  matrix  is  shown  in 
Table  A3,  while  the  C  matrix  is  in  Table  A4. 


TABLE  At 

THE  STATES  OF  THE  YF-17  STATE-SPACE  MODEL 


Z 

6 

4* -Structural  mode 
^-Structural  mode 
-Control  surface  mode 
42 -Control  surface  mode 
dZ/dt 
de/dt 
d4i/dt 
<U2/dt 
dd«/dt 


TABLE  A2 


THE  A  MATRIX  OF  THE  YF-iT 


ROW  1 

0 .OOOOE+OO  0 . OOOOE+OO  0. OOOOE+OO 
0.1000E+01  0. OOOOE+OO  0 . OOOOE+OO 
ROW  2 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0 . 1000E+01  0. OOOOE+OO 
ROW  3 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0.1000E+01 
ROW  4 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
ROW  S 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
ROW  6 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
ROW  7 

0. OOOOE+OO  0 . OOOOE+OO-O . 3640E+02 
-0. 1972E+01  0 .4510E+03-0 .3167E+00 
ROW  8 

0. OOOOE+OO  0. OOOOE+OO-O. 2S36E+00 
-0 . 5530E-0 1 -0 . 2834E+01 ~0 . 3462E-0 1 
ROW  9 

0. OOOOE+OO  0 .OOOOE+OO-O . 1063E+04 
-0 . 6238E+01  0. 1225E+01-0 .2530E+01 
ROW  10 

0. OOOOE+OO  0. OOOOE+OO-O. 2731E+03 
-0 . 9393E+01  0 .541 9E+ 01-0 . 1407E+01 
ROW  11 

0. OOOOE+OO  0. OOOOE+OO-O. 6130E+01 
-0 . 2628E+00-0 . 2952E+02-0 . 2426E+00 
ROW  12 

0. OOOOE+OO  0. OOOOE+OO  0.4492E+02 
0 . 2328E+00-0 . 3660E+02  0.3252E+00 


MODEL  AT  458  FT/SEC 


0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0 . 1000E+01  0. OOOOE+OO  0. OOOOE+OO 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0 . 1000E+01  0. OOOOE+OO 

0. OOOOE+OO  0. OOOOE+OO  0. OOOOE+OO 
0. OOOOE+OO  0. OOOOE+OO  0.1000E+01 

0.8012E+02  0 . 1733 E+ 03  0.6701E+03 
0.1976E-02  0 . 2958E+00  0.2355E+00 

0 .3157E+01-0 .5169E+02-0 .2357E+03 
0 . 6292E-0 1 -0 . 1 353E-01 -0 . 2464E+00 

0 . 5077E+03  0 . 8081E+03  0.5011E+04 
0 . 9873E+  00  0 . 6515E+00  0.3356E+01 

0 . 1428E+04  0 . 2240E+04  0.6060E+04 
0.2632E+01  0 . 2659E+0 1  0.2462E+01 

0 . 1664E+02-0 .9885E+05  0.2625E+02 
0 . 3485E+  00-0 .101 5E+0 1  0.U70E+00 

0 .4854E+02-0 . 1 162E+03-0 . 1395E+06 
0 .2543E+00-0 . 1 188E+00-0 . 2577E+02 


TABLE  A3 


THE  B  MATRIX  OF  THE  YF-17  MODEL  AT  458  FT/SEC 


0 .OOOOE+OO 
O.OOOOE+OO 
0. OOOOE+OO 
O.OOOOE+OO 
i  O.OOOOE+OO 
(  O.OOOOE+OO 
t  “0 . 2679E+03 
0. 1689E+03 
-0 . 1233E+04 
-0 . 2295E+04 
0 . 9885E+  05 
0.10UE+03 


0  .  OOOOE+OO" 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
O.OOOOE+OO 
-0.9205E+03 
0.3064E+03 
-0 .4638E+04 
-0.6524E+04 
-0 .3273E+02 
0 . 1396E+06 


TABLE  A4 

THE  C  MATRIX  OF  THE  YF-17  MODEL  AT  458  FT/SEC 
USING  TWIST  POSITION  SENSOR 


[0  0  0.008089  -0.01584  00000000] 


B. 


MAC  User's  Manual 


The  Methodology  for  Aeroelast ic ity  and  Controls  (MAC) 
program  was  developed  from  the  theory  presented  in  the  main 
body  of  this  thesis.  MAC  integrates  the  results  of  unsteady 
aerodynamics  codes,  structural  dynamics,  and  controls  into  a 
useful  tool  for  both  structural  dynamics  and  flight  control 
engineers.  MAC  was  developed  unt»er  AFIT  thesis  number 
AFIT/GAE/AA/86J-02  for  the  Flight  Dynamics  Laboratory.  The 
program  resides  on  the  Flight  Dynamics  Laboratory  Vax  11/785 
computer  with  a  VMS  operating  system.  MAC  accesses  the  IMSL 
library,  the  PLOTIO  plotting  library  and  the  DI3000  plotting 
library.  The  example  input  data  in  this  appendix  is  the 
YF-17  model  data  used  in  this  thesis. 

The  majority  of  input  to  MAC  can  either  be  input  from 
the  terminal  or  automatically  input  from  files.  There  are 
three  exceptions  to  this.  The  aerodynamic  forces  and  the 
print  input  data  must  be  input  from  files,  while  the 
feedback  matrices  must  be  input  from  the  terminal .  The 
following  logical  file  units  are  assigned  to  the  following 


data  files. 


TABLE  B1 


DATA  FILES  USED  IN  PROGRAM  MAC 

Logical  Unit  Data  File  Content 


1 

2 

3 

4 

5 

6 
T 

9 

10 
11 

12 

90 

98 

99 


Structural  Data 
Aerodynamic  Data 
Initial  Conditions 
Output  File 
Terminal  Input 
Terminal  Output 
Root  Locus  Input 
Output  Files  for  Use 
in  MATRIXx 
Print  Data 
Forcing  Function 
Matrix 
C  Matrix 

Data  for  Plotting 
Root  Locus 
Scratch  File 
Scratch  File 


The  first  thing  MAC  does  when  started  is  to  read  the 
system  print  parameters  from  a  file.  Then  the  user  will 
prompted  by  the  following  menu. 


be 


MAC  ready  for  input 
Enter  (#): 

1.  Automatic  input  of  items  (2)-<7) 

2.  Read  in  structural  data 

3 .  Read  in  aerodynamic  data 

4.  Read  in  initial  conditions 

5.  Read  in  root  locus  values 

6.  Read  in  forcing  function  matrix 

7.  Read  in  the  C  matrix 

8.  Completion  of  input  data 


If  the  user  chooses  option  i,  MAC  will  automatically  read  in 
all  the  data  from  items  2-7  from  files.  If  the  user  wishes 
to  input  the  data  separately,  the  options  2-8  must  be 
accomplished  in  order.  Once  the  data  is  input,  the 
remainder  of  the  program  will  prompt  the  user  for  input. 

The  terminal  input  is  in  free  field  format.  The  remaining 
portion  of  this  manual  describes  the  format  of  the  input 
data  on  files. 

Print  Data .  The  system  constants  for  the  determination 
of  what  is  to  be  output  is  in  this  file.  The  file  is  in 
namelist  format,  as  shown  in  Table  B2 .  Below  is  an 
explanation  of  each  variable  and  it’s  default  value. 

NAERO(O)  -  set  to  1  to  print  the  aerodynamic  matrices 
onto  unit  4,  the  output  file 


NMAT (0 ) 


set  to  1  to  print  the  X,  Y,  and  Z  matrices 
onto  unit  4 


NINV(O)  -  Set  to  1  to  print  the  X  1  matrix 
onto  unit  4 

NAMAT(O)  -  set  to  1  if  the  A  and  B  matrices  are  to  be 
set  to  files  that  are  compatible  with 
MATRIXx 

NEIG(l)  -  set  to  0  if  the  eigenvalues  are  not  to  be 
sent  to  the  plot  file 

NSTAB(O)  -  set  to  1  if  the  rigid  body  non-dimensional 
stability  derivatives  are  to  be  modified 


TABLE  B2 

PRINT  DATA  EXAMPLE 


$PRKT  NAEROrO,  HMAT= i,  NINV=0,  NAMAT= 1 ,  NEIG  =  1,  NSTAB = 0  $END 


Structural  Data .  The  structural  data  is  input  the 
diagonal  elements  of  mass  matrix,  and  the  associated 
frequency  (in  Hz)  of  vibration.  Table  B3  contains  an 
example.  The  data  is  input  in  the  following  format. 


Line  1  -  N=total  #  of  modes  (limited  to  20),  NRB= i  of 

rigid  modes  (limited  to  2),  N S=#  of  structural 
modes  (limited  to  13),  NC=#  of  control  surface 
modes  (limited  to  5) 

Format  (415) 

Line  2  -  Generalized  masses  in  the  order  of  rigid  modes, 
structural  modes,  and  control  surface  modes. 
Line  2  is  repeated  as  often  as  necessary  to 
account  for  all  the  modes 
Format  (4E15 . 7 ) 

Frequencies  of  the  modes  of  vibration  in  the 
same  order  as  above.  Line  3  is  also  repeated 
as  necessary 
Format  (4E15.7) 


Line  3 


TABLE  B3 


EXAMPLE  STRUCTURAL  DATA 


6  2  2 
0 .7061300E+01 
0 . 1027000E+00 
0 .OOOOOOOE+OO 
0 . 5000000E+02 


2 

0 .T776900E+02 
0. 1000000E+00 
0 .OOOOOOOE+OO 
0 . 6000000E+  00 


0 .3115000E+00 
0 .4628000E+01 


0 . 10 1 9000E+00 
0 .7186000E+01 


Aerodynamic  Data .  The  generalized  aerodynamic  forces 
(Q^jS),  from  unsteady  aerodynamic  codes  is  input  as  complex 
pairs  with  j  increasing  before  i.  Table  B4  is  an  example  of 
this  data  for  the  YF-17.  The  format  of  the  variables  is 
shown  below  Table  B4.  This  is  all  the  data  that  MAC  will  use 
from  this  file.  It  will  ignore  any  other  aerodynamic  data. 
Kl  should  be  0.0,  and  K2  should  be  a  small  value  of  reduced 
frequency . 


Initial  Condition  Data .  This  data  includes  the  density 
speed  of  sound,  reference  areas  and  lengths,  and  the  mode 
shape  magnitudes.  Table  B5  is  an  example  of  this  data  for 
the  YF-17.  The  variable  format  is  explained  below. 


3 

Line  1  -  RHOO=density  used  in  slugs/ft  ,  A0=speed  of 
sound  at  that  density  altltude^in  ft/sec, 

SREF -reference  wing  area  in  ft  , 

CBAR - ref erence  chord  (Mean  Aerodynamic  Chord) 
if  ft,  BR=reference  semi-chord  in  ft 
Format  (F10 .7,  Flo . 2,  3F10 . 5  ) 


Line  2  -  ALPHAO  and  HO,  the  magnitudes  of  the  pitch  and 
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si 


plunge  mode  shapes 
Format  (2F10.5) 

Line  3  -  ZETA(I)=the  magnitudes  of  the  structural  mode 
shapes.  This  line  is  repeated  as  often  as 
necessary  to  account  for  all  the  elastic  modes 
Format  (6F10.5) 

Line  4  -  DELTA( I ) = the  magnitudes  of  the  control  surface 
mode  shapes 
Format  (6F10.5) 


TABLE  B5 

EXAMPLE  INITIAL  CONDITION  DATA 


0.0010700  1100.00  13.68000  2.97000  1.48500 

12.00000  1.00000 
1.00000  1.00000 
12.00000  12.00000 


Root  Locus  Data .  The  root  locus  data  contains  the 
velocities  that  are  to  be  used  to  compute  the  eigenvalues 
for  plotting.  The  maximum  number  is  20.  Table  B6  shows  an 
example . 


Line  l  -  NVEL= »of  velocities  to  be  used 
Format  (15) 

Line  2  -  V(I)=each  velocity  in  ft/sec,  repeated  NVEL 
t  imes 

Format  (F10.3) 


TABLE  B6 


EXAMPLE  ROOT  LOCUS  DATA 


15 

0.000 
100.000 
200.000 
250.000 
300.000 
350 .000 
381 .000 
400.000 
420  .000 
440.000 
458.000 
470.000 
480.000 
490.000 
500.000 


Fore ing  Matrix.  The  forcing  matrix  is  used  to  create 
the  B  matrix.  It  is  defined  in  the  main  text  of  the 
thesis.  Table  B7  is  the  example  used  for  the  YF-17. 


Line  1  -  NINP=#of  inputs  the  system 
Format  (IS) 

Line  2  -  QMAT( I,  J ) =the  forcing  matrix  (NxNINP) 

input  row  first,  J  increases  before  i.  Line  2 
is  repeated  as  often  as  necessary. 

Format  (8E15.7) 


TABLE  BT 


EXAMPLE  P0RC1HG  MATRIX 


l 

I.NWNM8«M 

0.00000008*00 


0.0000000  E*00 
0 . 0000000E*00 


0. 808O000E-02 
0 . 00000008*00 


-0. 1684000E-01 

0.0000000E*00 


0 . 0000000E-00 
0.00000008*00 


0 . 0000000E*00 
0.00000008*00 


Measurement  Matrix .  The  C 
file.  Table  B8  is  the  example 


matrix  is  input  from  this 
used  for  the  YF-1T . 


Line  1  -  NOUT=#of  outputs  into  the  system 
Format  (15) 

Line  2  •?  CMAT ( I,  J )  =  the  measurement  matrix  (HOUTxH) 

input  rotr  first,  j  increases  before  1.  Line  2 
is  repeated  as  often  as  necessary. 

Format  (6ELS.T) 


TABLE  B8 

EXAMPLE  MEASUREMENT  MATRIX 

» 


2 

0 . 0000000B+00  0 . 0000000B+00  0.10000008*01  0.10000008*01  0.00000008*00  0.00000008*00 

0.00000008*00  0.00000008*00  0.10000008*01  0.00000008*00  0.00000008*00  0.10000008*01 


The  remain  inf  inputs  are  from  the  terminal,  as 
requested  by  the  the  program  prompts. 


This  Appendix  contains  the  program  listings  of  all  the 
programs  that  comprise  the  MAC  program.  MAC  uses  the  IMSL 
subroutine  library  for  matrix  inversion  (LINViF),  matrix 
multiplication  (VMULFF),  and  eigenvalue  computation  (EIGRF). 
MAC  also  uses  the  PLOTIO  plot  routine  library  in  the 
velocity  root  locus  plotting  program  to  clear  the  screen 
before  plotting.  Finally,  MAC  uses  the  DI3000  plot  routine 
library  to  create  both  screen  plots  and  the  hard  copy 
plots.  The  main  calling  program  (MAC)  is  listed  first.  The 
remaining  subroutine  libraries  are  in  alphabetical  order. 

All  the  routines  are  listed  below  as  they  appear  in  this 
append ix . 

TABLE  Cl 

MAC  PROGRAM  LISTING  ORDER 


MAC 

AMAT 

AUGMENT 

BMAT 

EVAL 

FASTCHG 


INPUT 

MATSAV 

MINV 

MMULT 

RLPLOT 

STABDER 


The  subroutines  are  fairly  well  commented  so  that  a  person 
familiar  with  the  theory  developed  in  the  main  text  and  with 


FORTRAN  should  be  able  to  interpret  the  program. 


PROGRAM  MAC 


C 

C  This  program  was  developed  under  AFIT  Thesis  AFIT/GAE/AA/86J-2 
C  by  John  J.  Cerra  II. 

C 

C  This  program  uses  generalized  mass  and  stiffness  matrices, 

C  along  with  generalized  forces  from  an  external  source  (FASTOP) 

C  and  combines  them  into  an  appropriate  state  space  formulation  for 
C  use  by  any  control  analysis  and  design  package.  The  state  space 
C  model  is  the  minimum  order  attainable,  2N,  where  N  is  the  number 
C  of  modes  used  (including  rigid  body  and  control  surface  modes). 

C 

C  Subroutine  INPUT  reads  in  the  generalized  mass  and  stiffness 
C  matrices.  It  also  reads  in  the  generalized  forces  from  FASTOP  at 
C  two  reduced  frequencies  (usually  k=0  and  a  small  k) .  The  initial 
C  conditions  are  also  input  here,  along  with  the  root  locus  data. 

C  Finally  it  also  inputs  the  forcing  function  mat.!x,  Q, 

C  and  the  C  matrix  if  needed. 

C 

C  Subroutine  STABDER  computes  stability  derivatives  from  the  FASTOP 
C  unsteady  aerodynamics.  There  is  also  an  option  to  modify  the 
C  derivatives  if  derivatives  from  a  better  source  are  known. 

C 

C  Subroutine  AMAT  uses  the  above  data  to  compute  the  state  space 
C  A  matrix  used  for  stability  analysis.  The  order  of  A  is  2N.  A 
C  is  left  as  a  function  of  dynamic  pressure  so  that  a  velocity 
C  root  locus  can  be  performed  as  in  a  typical  flutter  solution. 

C 

C  Subroutine  BMAT  computes  the  B  matrix  from  the  forcing  function 
C  matrix  Q.  B  is  a  function  of  dynamic  pressure. 

C 

C  Subroutine  MINV,  called  by  AMAT,  computes  the  inverse  of  the  X 
C  matrix  using  an  IMSL  routine  LINV1F. 

C 

C  Subroutine  MMULT,  called  by  AMAT,  and  BMAT  multiplies  matrices. 

C  This  is  used  to  form  the  A  and  8  matrices. 

C 

C  Subroutine  EVAL  computes  the  eigenvalues  of  the  A  matrix 
C  The  eigenvalues  for  various  velocities  can 

C  be  accomplished  and  then  can  be  plotted  on  the  real  and  imaginary 
C  plane  to  produce  a  velocity  root  locus. 

C 

C  Subroutine  AUGMENT  augments  the  A  matrix  with  the  appropriate 
C  feedback,  to  form  the  augmented  (closed  loop)  A  matrix 
C 

C  Subroutine  FASTCHG  creates  non-dimensional  coefficients  out  of 
C  the  unsteady  aerodynamics  of  FASTOP 
C 

C  Subroutine  RLPLOT  plots  the  eigenvalues  of  the  A  matrix  with  respect 
C  to  velocity  on  the  complex  plane,  creating  a  ’velocity  root  locus’. 


Subroutine  MATSAV  saves  matrices  to  a  file  in  a  compatible 
form  for  MATRIXx 

CHARACTER *1  QAUG , QLPT 
CHARACTER*80  TITLE 
INTEGER  NVEL 
REAL  VE,V(20) 

C0MM0N/RL/NVEL,V 

Open  output  file 

OPEN (4 , STATUS*  * NEW  * ) 

CALL  INPUT 

OPEN (98 , STATUS* » SCRATCH  * ) 

WRITE (98,*)  NVEL 

Velocity  root  locus  loop  for  the  unaugmented  aircraft 

IF  (NVEL.GT.l)  THEN 
DO  1  1=1, NVEL 
CALL  STABDER (NVEL, V (I) ) 

CALL  AMAT 
CALL  BMAT (V (I) ) 

CALL  EVAL 
1  CONTINUE 
ELSE 
NVEL=0 
VE=0 . 001 

CALL  STABDER  (NVEL,  VE) 

CALL  AMAT 
CALL  BMAT (VE) 

CALL  EVAL 
ENDIF 

Option  to  plot  the  open  loop  aircraft  roots 

WRITE(*, ’ (A) ’)  *  Do  you  want  to  plot  the  open  loop  roots?’ 
READ (*,100)  QPLT 
IF  (QPLT.EQ. ’Y’)  THEN 

WRITE(*, ’ (A) ’)  ’  What  title  do  you  want  for  the  plot?’ 
READ (*,100)  TITLE 

CALL  RLPLOT (0.0, 8. 0,-6. 0,6.0, TITLE) 

CLOSE (90) 

ENDIF 

WRITE(* , ’ (A) ')  ’  Do  you  want  to  augment  the  A  matrix?’ 

READ (*,100)  QAUG 
IF  (QAUG.EQ. ’Y’)  THEN 

Call  to  closed  lop  augmentation  program 


§ 


CALL  AUGMENT 
ENDIF 
IDAT=99 
CLOSE(l) 

CLOSE (2) 

CLOSE (3) 

CLOSE (4) 

CLOSE (7) 

CLOSE (90) 

CLOSE (9) 

CLOSE (10) 

CLOSE(ll) 

CLOSE (12) 

CLOSE (98) 

WRITE(* , ’ (A) ’)  ’  This  program  has  ended 
100  FORMAT  (A) 

END 
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SUBROUTINE  AMAT 
C 

C  This  subroutine  uses  the  stability  derivatives  above  for  use  in 
C  developing  the  open  loop  reduced  order  A  matrix  for  stability 
C  analysis,  and  control  work.  It  is  of  much  smaller  order  than 
C  that  resulting  from  using  Pade  fit  aerodynamics. 

C 

CHARACTER* 1  QASAV 

CHARACTER *20  ANAM 

INTEGER  N , NRB , NS , NC , NV 

COMMON / INDEX/N , NRB , NS , NC 

INTEGER  NAERO , NMAT , NINV , NAMAT , NEIG 

COMMON/PRT /NAERO , NMAT , NINV , NAMAT , NEIG 

REAL  ZALPHA , MALPHA , ZALPHDT , MALPHDT ,ZQ,MQ, 

OFZALPHA , FZALPDTN , FZETQ , ZZETA , MZETA , 

OZZETAOT, MZETAOT, ZZETDDT, MZETDDT, 

OFZETZET , FZETZDT , FZEZDDT , 

OFZETDEL , FZETDDT , FZEDDDT , 

OFDELZET , FDELZDT , FDEZOOT , 

OZDELTA ,  MDaTA ,  FDALPHA ,  FDALPOT ,  FDETQ , 

OZDaTDT ,  MDELTDT ,  ZDaDDT ,  MDELDDT , 

OFOaOEL ,  FDELDDT,  FDEDOOT,  VEL 

DIMENSION  FZALPHA(15) ,FZALPDT(15) ,FZETQ(15) ,ZZETA(15) ,MZETA(15) , 
OZZETADT (15) ,MZETADT (15) , ZZETDOT (15) .MZETDDT (15) , 

OFZETZET (15,15), FZETZDT (15,15), FZEZDDT (15 , 15) , 

OFZETDEL (15, 5) , FZETDDT (15, 5) , FZEDDDT (15, 5) , 

OFDELZET (5,15), FDELZDT (5, 15) ,FDEZDDT (5,15) , 

0ZDaTA(5)  ,  MDELTA  (5)  ,  FDALPHA  (5)  ,  FDALPDT  (5)  ,FDETQ(5)  , 

OZDaTDT  (5)  ,  MDaTDT  (5)  ,  ZDaDDT  (5)  ,  MDELDDT  (5)  , 

OFDaDEL  (5,5),  FDELDOT  (5,5),  FDEDDDT  (5,5) 

COMMON/DERI V/ZALPHA , MALPHA , ZALPHDT , MALPHDT ,ZQ,MQ, FZALPHA , FZALPDT , 
OFZETQ , ZZETA , MZETA , ZZETADT , MZETADT , ZZETDDT , MZETDDT , FZETZET , 
OFZETZDT,  FZEZDDT,  FZETDa,  FZETDDT,  FZEDDDT,  FDELZET,  FDELZDT,  FDEZDDT, 
OZDaTA ,  MDaTA ,  FDALPHA ,  FDALPDT ,  FDETQ , 

OZDaTDT ,  MDaTDT ,  ZDaDDT ,  MDELDDT ,  FDELDEL ,  FDELDDT ,  FDEDDDT ,  VEL 
REAL  GM,GK 

DIMENSION  GM(20) ,GK(20) 

COMMON/STRUCT /GM , GK 
REAL  K1,K2 
COMPLEX  GF1,GF2 

DIMENSION  GF1(20,20),GF2(20,20 
COMMON/AERO /K1 , K2 , GFl , GF2 , PI 
REAL  X , Y , Z , AMATRIX , XINV , XINVY , XINVZ 

DIMENSION  X (20 , 20) , Y (20 , 20) , Z (20 , 20) , AMATRIX (40 , 40) , XINV (40) , 
OXINVY (20 , 20) , XINVZ (20 , 20) 

COMMON/INV/X,XINV 
COMMON/MMUL/Y , Z , XINVY , XINVZ 
COMMON/EIG/AMATRIX 
REAL  XM0M(5) 

COMMON/MOMENT/XMOM 


t 

1 

v 

V 

V 

V 

V 
‘4 
‘4 


V, 

A 

i 

i 


3 
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■>.  y  y.  y  »•.  vv.v-.v 1 


000  ooo  ooo  ooo  ooo 


00  2  1=1, NS 
J=I+NRB 

IF  (NRB.EQ.O)  GOTO  20 
X(J, 1)=-FZALPDT (I) /VEL 
X (1 , J)=-ZZETDDT (I) 

X(2, J)=-MZETDDT (I) 

Y (J, 1)=-FZALPHA(I) /VEL 
Y ( J / 2) =-FZETQ (I) 

Y (1 , J)=-ZZETADT (I) 

Y (2 , J) s-MZETAOT (I) 

Z(l, J)=-ZZETA(I) 

Z(2, J)=-MZETA(I) 

20  CONTINUE 
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IF  (NS.EQ.O)  GOTO  4 

Elastic/Elastic  aero 

DO  4  K=1 , NS 
L=K+NRB 

X ( J , L) =-FZEZDDT (I,K)+X(J,L) 

Y ( J , L) =-FZETZDT (I , K) 

Z ( J , L) =-FZETZET (I , K) +Z ( J , L) 

4  CONTINUE 
IF  (NC.EQ.O)  GOTO  21 

Elastic/Control  and  Control/Elastic  aero 

DO  5  K=1 , NC 
L=K+NS+2 


VZVUW  V"  JTTJTV 


X  (L ,  J) =-FDEZDDT (K , I) 
X(J,L)=-FZEDDDT (I,K) 
Y(L,J)=-FDELZDT(K,I) 
Y(J,L)=-FZETDDT(I,K) 
Z  (L ,  J) =-FDELZET (K , I) 
Z  ( J ,  L) =-FZETDEL (I , K) 
5  CONTINUE 
21  CONTINUE 
2  CONTINUE 

IF  (NRB.EQ.O)  GOTO  23 
IF  (NC.EQ.O)  GOTO  22 


'.VI 


M 


v' 
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Control /Rigid,  Rigid/Control  aero 

DO  3  1=1 ,NC 
J=I+NRB+NS 

X(J, 1)=-FDALPDT (I)/VEL 
X (1 , J) =-ZDELDDT (I) 

X (2 , J) =-MDELDDT (I) 

Y (J, 1)=-FDALPHA(I) /VEL 
Y(J,2)=-FDETQ(I) 

Y (1 , J) s-ZDELTDT (I) 

Y (2 , J) =-MDELTDT (I) 

Z (1 , J)=-ZDELTA(I) 

Z(2, J)=-MDELTA(I) 

23  CONTINUE 

IF  (NC.EQ.O)  GOTO  6 
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Control /Control  aero 


DO  6  K=1 ,NC 
L=K+NS+NRB 

X ( J , L) =-FDEDDDT (I,K)+X(J,L) 
Y(J,L)=-FDELDDT(I,K) 

Z ( J , L) =-FDELDEL (I , K) +Z  ( J , L) 
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6  CONTINUE 
3  CONTINUE 
DO  32  J=1,NC 
L=J+NS+NRB 
XM0M(J)=Z(L,L) 

32  CONTINUE 
22  CONTINUE 

IF  (NMAT.EQ.O)  GOTO  14 

WRITE (4, ’(/A) ’)  ’  The  X  matrix  is’ 

DO  10  1=1, N 
WRITE (4, 300)  I 
DO  11  J=1,N,6 

WRITE (4, 301)  (X(I,K) ,K=J, J+5) 

11  CONTINUE 
10  CONTINUE 

WRITE(4,  ’(/AV)  >  The  Y  matrix  is’ 

DO  12  1=1, N 
WRITE (4, 300)  I 
DO  13  J=1,N,6 

WRITE (4, 301)  (Y (I,K) ,K=J, J+5) 

13  CONTINUE 

12  CONTINUE 

WRITE (4, ’(/A) ’)  ’  The  Z  matrix  is’ 

DO  14  1=1, N 
WRITE (4, 300)  I 
DO  15  J=1 ,N,6 

WRITE (4,301)  (Z (I , K) ,K=J, J+5) 

15  CONTINUE 

14  CONTINUE 

Initialization  of  the  A  matrix. 

DO  31  1=1,  N2 
DO  31  J=1,N2 
AMATRIX (I , J)  =0 . 0 
31  CONTINUE 

Call  to  invert  X  and  to  multiply  XINV*Y  and  XINV*Z 

CALL  MINV 
CALL  MMULT 
DO  8  1=1, N 
DO  8  J=1,N 
K=I+N 
L=J+N 

AMATRIX (K , L) =-XINVY (I , J) 

AMATRIX ( I, K) =1.0 
AMATRIX (K,J)=-XINVZ (I, J) 

8  CONTINUE 
IF  (NAMAT.EQ. 1)  THEN 
OPEN (99 , STATUS= ’ SCRATCH ’ ) 


nnn 


NV=NINT (VEL) 

WRITE (99 ,310)  NV 
310  FORMAT ( ’ [ . MATRXX] A ’ , 14 . 4) 
REWINO  99 
READ (99, 302)  ANAM 
CLOSE (99) 

DO  40  1=1, N2 
DO  40  J=1,N2 
AI J=AMATRIX (I , J) 

A(I,  J)=DBLE(AIJ) 

40  CONTINUE 

OPEN (9 , FILE=ANAM , STATUS= ’ NEW ’ ) 


Call  to  the  routine  to  put  the  A  matrix  in  MATRIXx  format 


CALL  MATSAV (9 , ANAM , 40 , N2 , N2 , 0 , A , DUMMY , ’ (1P8E15.7) ’) 
CLOSE (9) 

ENDIF 

WRITE (4, ’ (/A) ’)  ’  The  A  matrix  of  dx/dt=A*x  is’ 

DO  16  1=1, N2 
WRITE (4, 300)  I 
DO  17  J=1 ,N2,6 

WRITE (4 , 301 )  (AMATRIX (I,K) ,K=J, J+5) 

17  CONTINUE 
16  CONTINUE 

WRITE (98, *)  ( (AMATRIX (I , J) , J=1 , N2) , 1=1 , N2) 

IF  (NRB.NE.O)  WRITE (4, 304) 

IF  (NS.NE.O)  WRITE (4, 305)  (I, 1=1, NS) 

IF  (NC.NE.O)  WRITE (4, 306)  (I,I=1,NC) 

IF  (NRB.NE.O)  WRITE (4, 307) 

IF  (NS.NE.O)  WRITE (4, 308)  (I, 1=1, NS) 

IF  (NC.NE.O)  WRITE (4, 309)  (I,I=1,NC) 

300  FORMAT ( *  ROW’, 13) 

301  FORMAT (IX, 6E1 1.4) 

302  FORMAT (A) 

304  FQRMAT(/,’  The  states  of  the  A  matrix  are:’,/, 

O’  Z’  / 

O’  THETA’) 

305  FORMAT (’  ZETA( ’ , 12, ’) -Structura I  mode’) 

306  FORMAT(’  DELTA( ’, 12, ’) -Control  surface  mode’) 

307  FORMAT (’  dZ/dt’,/,’  dTHETA/dt’) 

308  FORMAT ( ’  dZETA(’ ,12, ’)/dt’) 

309  FORMAT (’  dDELTA(’,I2, ’)/dt’) 

303  FORMAT (6E15. 7) 

RETURN 
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SUBROUTINE  AUGMENT 

This  subroutine  will  augment  the  state  space  equations  with 
the  appropriate  feedback  (output  or  state)  to  form  the  closed 
loop  state  equations 

REAL  AUG (40 , 40) , ACL (40 , 40) , TEMP (40 , 40) , CMAT (40 , 40) , 
OBMATRIX (40 , 40) , FMAT (5 , 40) , V (20) 

INTEGER  N,N2,NRB,NS,NC, IFDBK , NP , NQUT , I ER , NA , NB 
CHARACTER *1  QPLT,QCHG 
CHARACTER* 80  TITLE 
COMMON/EIG/ACL 
COMMON/INDEX/N , NRB , NS , NC 
COMMON/CMATRX/NOUT, CMAT 
COMMON/BMA  rX/NP , BMATRIX 
C  OMMON / RL / N VEL , V 
N2=2*N 
3  CONTINUE 
DO  8  1=1,5 
DO  8  J=1 , 40 
FMAT (I, J) =0.0 
8  CONTINUE 

Find  out  whether  output  or  state  feedback  is  to  be  used 

WRITE (* , 700) 

READ (*,701)  IFDBK 
IF  (IFDBK. GT. 2)  GOTO  3 
IF  (IFDBK. EQ.O)  GOTO  7 
IF  (IFDBK. Eq.l)  THEN 

output  feedback  augmentation  loop,  creates  BKC=AUG 

WRITE (*,702)  NP , NOUT 
NA=NP 
NB=NOUT 
DO  1  1=1, NP 
WRITE (*,703)  I 

READ(*,*)  (FMAT (I , J) , J=1 , NOUT) 

1  CONTINUE 
GOTO  13 
14  CONTINUE 

Call  to  IMSL  routine  to  multiply  B*K,  and  then  BK*C 

CALL  VMULFF (BMATRIX , FMAT , N2 , NP , NOUT  ,40,5, TEMP , 40 , 1 ER) 
CALL  VMULFF (TEMP , CMAT , N2 , NOUT , N2 , 40 , 40 , AUG , 40 , I ER) 

WRITE (4,  ’(A)  ’)  ’  For  OUTPUT  feedback’ 

WRITE(4, ’ (A)  ’)  ’  With  a  feedback  matrix  of’ 

DO  16  1=1, NP 
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WRITE (4,703)  I 
00  17  J=1,N0UT,6 
WRITE (4 ,706)  (FMAT (I ,  K) , K= J , J+5) 
17  CONTINUE 
16  CONTINUE 
ELSE 


state  feedback  augmentation  loop,  creates  BK=AUG 

WRITE (* , 702)  NP,N2 
NA=NP 
NB=N2 

DO  2  1=1, NP 
WRITE (*,703)  I 
READ(*, *)  (FMAT (I , J) , J=1 , N2) 

2  CONTINUE 
GOTO  13 
15  CONTINUE 

Call  to  IMSL  routine  to  multiply  B*K 

CALL  VMULFF (BMATRIX , FMAT , N2 , NP , N2 , 40 , 5 , AUG , 40 , IER) 
WRITE (4, ’(A) ’)  ’  For  STATE  feedback’ 

WRITE (4, ’ (A) ’)  ’  With  a  feedback  matrix  of’ 

DO  18  1=1, NP 
WRITE (4, 703)  I 
DO  19  J=1,N2,6 

WRITE (4, 706)  (FMAT (I , K) , K=J , J+5) 

19  CONTINUE 
18  CONTINUE 
ENDIF 

DO  4  J=1,N2 
DO  4  1=1, N2 
TEMP (I, J) =0.0 
4  CONTINUE 

C  Unit  98  contains  the  unaugmented  A  matrices  created  duri 
C  the  first  root  locus.  This  section  now  creates  another 
C  root  locus  loop,  computing  eigenvalues  and  using  the 
C  c I osed  I oop  A  matr i x  Ac  I =A-AUG 
C 

REWIND  (98) 

READ (98,*)  NVEL 
DO  5  K=1 , NVEL 

READ (98 , *)  ( (TEMP (I , J) , J=1 , N2) , 1=1 , N2) 

DO  6  1=1, N2 
DO  6  J=1 ,N2 

ACL (I , J) =TEMP (I , J) -AUG (I , J) 

6  CONTINUE 

WRITE (*,704)  V(K) 

WRITE (4, 704)  V(K) 
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CALL  EVAL 
5  CONTINUE 


Option  to  plot  the  root  locus  points 

WRITE(*, ’ (A) ’)  ’  Do  you  want  to  plot  the  closed  loop  roots?’ 
READ (*,705)  QPLT 
IF (QPLT. EQ. ’Y’)  THEN 

WRITE(* , ’ (A) ’)  ’  What  title  do  you  want  for  the  plot?’ 

READ (*,705)  TITLE 

CALL  RLPLOT  (0.0, 8. 0,-6. 0,6.0, TITLE) 

CLOSE (90) 

ENDIF 
GOTO  3 

This  section  gives  the  user  the  option  to  change  the  values 
of  any  portion  of  the  feedback  matrix. 

13  CONTINUE 
DO  10  1=1, NA 
WRITE (*,703)  I 
DO  11  J=1 , NB , 6 

WRITE (*,706)  (FMAT (I , K) , K= J , J+5) 

11  CONTINUE 
10  CONTINUE 

WRITE(*, ’ (A) ’)  ’  Do  you  want  to  change  any  values?' 

READ (*,705)  QCHG 
IF  (QCHG.EQ. ’N’)  GOTO  12 
9  CONTINUE 

WRITE(*, ’ (A) ’)  ’  Please  input  the  row,  column,  and  new  value.’ 
READ(*, *)  I, J, FMAT (I, J) 

WRITE(* , ’ (A) ’)  ’  Do  you  want  to  change  another  value?’ 

READ (*,705)  QCHG 
IF  (QCGG.EQ. ’Y’)  GOTO  9 

12  CONTINUE 

IF  (IFDBK.EQ.l)  GOTO  14 
IF  (IFDBK.EQ.2)  GOTO  15 

700  FORMAT (’  What  type  of  feedback  are  you  using?’,//, 


0’ 

0. 

Return  to  main  program’,/. 

0’ 

1. 

Output’ ,/, 

0’ 

2. 

State’ ,/) 

701  FORMAT (15) 

702  FORMAT(’  The  feedback  matrix  will  have’, 13,’  rows  by’, 

013,’  columns.’,/,’  Please  input  the  feedback  matrix  by  rows.’,/) 

703  FORMAT (’  Row ’,13) 

704  FORMAT (//,’  The  closed  loop  eigenvalues  for  ’,F10.3,’  ft/sec’) 

705  FORMAT (A) 

706  FORMAT (IX, 6E1 1.4) 

7  CONTINUE 

RETURN 

END 
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SUBROUTINE  BMAT (VEL) 


This  subroutine  takes  the  forcing  function  matrix  and  converts  it 
into  the  B  matrix  of  dx/dt=Ax+Bu  where  u  is  the  control  input  vector 

REAL  X (20 , 20) , XINV (20 , 20) , QMAT (20,5), BMATRIX (40 , 5) , BMATl (20 , 5) 
INTEGER  N , NRB , NS , NC , NINP , NAERO , NMAT , NINV , NAMAT , NEIG , N2 , NV 
REAL*8  B (40, 40), DUMMY 
CHARACTER *20  BNAM 

COMMON/PRT /NAERO , NMAT , NINV , NAMAT , NEIG , NV 

COMMON /INDEX /N , NRB , NS , NC 

COMMON/BMATRX/NINP , QMAT 

COMMON/INV/X,XINV 

COMMON/BMATX/NP , BMATRIX 

REAL  XM0M(5) 

COMMON/MOMENT/XMOM 

Call  to  IMSL  routine  to  multiply  XINV*QMAT 

CALL  VMULFF (XINV , QMAT ,N,N, NINP ,20,20, BMATl , 20 , IER) 

NP=NINP 

Creates  the  actual  B  matrix 

DO  2  1=1, N 
K=N+I 

DO  2  J=1 ,NINP 

BMATRIX (K , J) =BMAT1 (I , J) *XMOM ( J) 

2  CONTINUE 
WRITE (4, 603) 

N2=2*N 
DO  1  1=1, N2 
WRITE (4, 600)  I 

WRITE (4 , 601)  (BMATRIX (I, J) , J=l, NINP) 

1  CONTINUE 
IF  (NAMAT. EQ.l)  THEN 
OPEN (99 , STATUS= ’ SCRATCH ’ ) 

NV=NINT (VEL) 

WRITE (99, 604)  NV 

604  FORMAT ( ’ [ . MATRXX] B ’ , 14 . 4) 

REWIND  99 

READ (99, 605)  BNAM 

605  FORMAT (A) 

CLOSE (99) 

DO  10  1=1, N2 
DO  10  J=1,NINP 
BI J=BMATRIX (I , J) 

B(I, J)=DBLE(BIJ) 

10  CONTINUE 
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OPEN (9 , FILE=BNAM , STATUS= ’ NEW  * ) 


Cali  to  routine  that  puts  the  6  matrix  into  MATRIXx  format 


CALL  MATSAV(9,BNAM,40,N2,NINP,0,B,DUMMY,  ’  (1P8E15.7)  ») 
CLOSE (9) 

ENOIF 

600  FORMAT (’  Row’, 13) 

601  FORMAT (1X,6E12.4) 

602  FORMAT (6E15. 7) 

603  FORMAT (//,’  The  B  matrix  of  dx/dt=Ax+Bu  is’) 

RETURN 

END 


-74- 


SK* 


PS 

i 


iLTi.t 

HLa’I 


I 


& 


$ 

1 


I 

if 


r>  r>  o 


SUBROUTINE  EVAL 


C  This  subroutine  computes  the  eigenvalues  of  the  A  matrix 
C  using  the  IMSL  routine  EIGRF. 

C 

REAL  AMATRIX,WK 

DIMENSION  AMATRIX (40 , 40) , WK (40) 

COMMON/EIG/AMATRIX 
COMPLEX  EIGVAL(40) 

INTEGER  N , NRB , NS , NC , IER , N2 
COMMON/INDEX/N , NRB , NS , NC 
INTEGER  NAERO , NMAT , NINV , NAMAT , NEIG 
COMMON/PRT /NAERO , NMAT, NINV , NAMAT, NEIG 
CHARACTER* 1  QEVAL 
CHARACTER *64  EIGFILE 
REAL*8  XX, XY ,XZ 
N2=N*2 
PI=3. 14159 

Call  to  IMSL  eigenvalue  solver  routine 

CALL  EIGRF (AMATRIX , N2 , 40 , 0 , EIGVAL , Z , IZ , WK , IER) 

IF  (IER.NE.O)  THEN 
WRITE (*,600)  IER 
STOP 
ENDIF 

IF  (NEIG.EQ.l)  OPEN (90, STATUS*’ SCRATCH’, FORM=’ UNFORMATTED’) 
WRITE(*. ’(/A) *)  »  The  eigenvalues  are:’ 

WRITE (*, ’  (A)  ’)  '  I  Eigenvalue  Freq(Hz) ’ 

WRITE(4, ’ (/A) ’)  ’  The  eigenvalues  are:’ 

WRITE (4, ’  (A)  ’)  »  I  Eigenvalue  Freq(Hz)’ 

DO  2  1=1, N2 

FREQ=AIMAG (EIGVAL (I) ) / (2 . *PI) 

XY=DBLE(FREQ) 

RV=REAL (EIGVAL (I)) 

XX=DBLE (RV) 

XIM=AIMAG (EIGVAL (I)) 

XZ=DBLE (XIM) 

WRITE (* , 602)  I , REAL (EIGVAL (I) ) , AIMAG (EIGVAL (I) ) , FREQ 
WRITE (4 , 602)  I , REAL (EIGVAL (I) ) , AIMAG (EIGVAL (I) ) , FREQ 
IF  (NEIG.EQ.l)  WRITE (90)  XX,XY,XZ 
2  CONTINUE 
1  CONTINUE 
WRITE (4, ’(A) ’)  ’1’ 

600  FORMAT(///, *****  There  was  an  error  in  the  eigenvalue’, 

O’  computation,  f  ’,13,’  ****’,///) 

601  FORMAT (A) 

602  FORMAT (I3,2Ell .4, ’ i ’ ,E11 .4) 

RETURN 
END 


SUBROUTINE  FASTCHG 


ft 


C  This  subroutine  changers  FASTOP  generalized  forces  into 
C  non-dimensional  form  by  multiplying  by  s**2  and  by  dividing 
C  by  SREF  for  forces  and  by  SREF*CBAR  for  moments 
C  also  the  difference  in  the  coordinate  systems  is  taken  into 
C  account.  FASTOP  is  a  left-hand  rule  axis  system  with  positive 
C  Z  up  and  positive  moment  down.  Normal  stability  in  body  axes 
C  positive  Z  is  downward  and  positive  moment  is  up. 

C 

COMPLEX  GF1,GF2 

REAL  S , SREF , CBAR , PI , SR , CB , Kl , K2 , VO , RHOO , ALPHAO , HO , ZETAO, 
ODELTAO.BR 
INTEGER  N,NRB,NS,NC 

DIMENSION  GF1 (20,20) ,GF2(20, 20) , ZETA0(15) ,DELTA0(5) 
C0MM0N/AER0/K1 , K2 , GF1 , GF2 , PI 

COMMON/IC /VO , RHOO , SREF , CBAR , ALPHAO , HO , ZETAO , DELTAO , BR 

COMMON/ INDEX/N , NRB , NS , NC 

NRB1=NRB*1 

NSNRB=NS+NRB 

NSNRB1=NSNRB*1 

NCNSNRB=NSNRB+NC 

S=12.0 

SR=SREF* 144.0 
CB=CBAR*12.0 
IF  (NRB.EQ.O)  GOTO  1 
DO  1  J=1,N 

GF1  (1 , J)=-GFl (1 , J) *S**2/SR 
GF2 (1 , J) =-GF2 (1 , J) *S**2/SR 
GF1 (2, J) =-GFl (2 , J) *S**3/ (SR*CB) 

GF2(2, J)=-GF2(2, J) *S**3/(SR*CB) 

1  CONTINUE 

IF  (NS.EQ.O)  GOTO  2 
DO  2  I=NRB1 ,NSNRB 
DO  2  J=1,N 

GFl (I, J)=-GF1 (I, J) *S**2/SR 
GF2(I , J)=-GF2(I , J) *S*+2/SR 

2  CONTINUE 

IF  (NC.EQ.O)  GOTO  3 
DO  3  I=NSNRB1 , NCNSNRB 
DO  3  J=1,N 

GFl (I , J) =-GFl (I , J) *S**3/ (SR^CB) 

GF2(I , J)=-GF2(I , J) *S**3/ (SR*CB) 

3  CONTINUE 
RETURN 
END 


* 
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SUBROUTINE  INPUT 


This  subroutine  contains  all  the  input  necsssary  to  run  the  program. 
The  total  dimensions  are  20  total  modes,  2  rigid  body  modes,  5  control 
surface  modes  and  15  elastic  modes  maximum. 


INTEGER  N , NRB , NS , NC , NSQ , NVEL , INPT , NINP , N2 , NOUT 

INTEGER  NAERO , NMAT , NINV , NAMAT , NEIG , NSTAB 

CHARACTER *64  STRUCT , AERO , INCOND , RTLCSF , INPTM .TITLE , OUTPM 

CHARACTER *1  Ql , Q2 , Q3, Q4 , Q5 , Q6 , Q7 , Q8, Q9 , QlO, Ql 1 

REAL  K 1 , K2 , VO , RHOO , SREF , CBAR , ALPHAO , HO , XMACH , AO , BR 

REAL  GM , GK , ZETAO, DELTAO , OMEGA , V , F (20 , 20) , G (20 , 20) 

REAL  QMAT (20,5), CMAT (40,40) 

DIMENSION  GM(20) , GK (20) , ZETAO (15) .DELTAO (5) .OMEGA (20) ,V(20) 
COMPLEX  GF1.GF2 

DIMENSION  GFl (20,20) ,GF2(20,20) 

COMMON/STRUCT /GM , GK 
C0MM0N/AER0/K1 , K2 , GFl , GF2 , PI 
COMMON/INDEX/N , NRB , NS , NC 

COMMON/ IC/ VO , RHOO , SREF , CBAR , ALPHAO , HO , ZETAO , DELTAO , BR 
COMMON/RL/NVEL , V 

C  OMMON/PRT /NAERO , NMAT , NINV , NAMAT , NE I G 
COMMON/BMATRX/NINP , QMAT 
COMMON/CMATRX /NOUT , CMAT 

NAMELIST /PRNT /NAERO , NMAT , NINV , NAMAT , NEIG , NSTAB 


Default  values  for  the  namelist  variables 

NAER0=O 

NMAT-0 

NINV=0 

NAMAT=0 

NEIG=1 

NSTAB=0 

OPEN(10,STATUS=’OLD’) 

READ (10, PRNT) 

PI=3. 14159 


Question  to  find  out  whether  the  data  is  to  be  automatically 
input,  or  whether  it  will  be  done  manually. 

24  CONTINUE 
WRITE (*,120) 

READ (*,121)  INPT 
IF  (INPT.EQ. 1)  THEN 

Auto  input  of  Structural  info 

OPEN (1,STATUS=* OLD’) 

READ(1 ,116)  N,NRB,NS,NC 
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N2=2*N 

READ(1,117)  (GM(I) , 1=1 ,N) 
READ (1,117)  (OMEGA (I), 1=1 ,N) 


Auto  i nput  of  Aero 

OPEN (2 , ST ATUS= ’ OLD ’ ) 

READ (2, 111)  TITLE 
READ (2, 128)  N , XMACH 
READ (2 ,127)  K1 

READ(2, 126)  (((F(I, J) ,G(I, J)) ,I=1,N) , J=1,N) 
DO  7  1=1, N 
DO  7  J=1,N 

GF1 (I , J) =CMPLX (F (I , J) , G (I , J) ) 

7  CONTINUE 

READ (2, 127)  K2 

READ (2, 126)  (((F(I, J) ,G(I, J)) ,I=1,N) , J=1,N) 
DO  11  1=1, N 
DO  11  J=1 ,N 

GF2 (I , J) =CMPLX (F (I , J) , G (I , J) ) 

11  CONTINUE 


Auto  input  of  initial  conditions 

OPEN (3 , STATUS= ’ OLD ’ ) 

READ (3, 118)  RHOO , AO , SREF , CBAR , BR 
IF  (NRB.NE.O)  READ(3, 119)  ALPHAO, HO 
IF  (NS.NE.O)  READ (3, 119)  (ZETAO(I) ,1=1, NS) 
IF  (NC.NE.O)  READ(3, 119)  (DELTAO(I) ,I=1,NC) 


Auto  input  of  root  locus  info 


OPEN (7 , STATUS* ’ OLD ’ ) 
READ (7, 121)  NVEL 
DO  15  1=1, NVEL 
READ (7, 122)  V(I) 

15  CONTINUE 


Auto  input  of  forcing  function  matrix 


OPEN (11,STATUS=’ OLD’) 

READ(11,121)  NINP 

READ (11, 123)  ((QMAT(I,J) ,J=1, NINP), 1=1, N) 


Auto  input  of  C  matrix 


OPEN ( 1 2 , ST ATUS= ’ OLD ’ ) 

READ (12, 121)  NOUT 

READ (12 , 123)  ( (CMAT (I , J) , J=1 , N2) , 1=1 , NOUT) 
GOTO  23 
ENDIF 


IF  (INPT.LE.O)  GOTO  24 
IF  (INPT.EQ.2)  GOTO  20 
IF  (INPT.EQ.3)  GOTO  21 
IF  (INPT.EQ.4)  GOTO  22 
IF  (INPT.EQ.5)  GOTO  26 
IF  (INPT.EQ.6)  GOTO  25 
IF  (INPT.EQ.7)  GOTO  27 
IF  (INPT.EQ.8)  GOTO  23 
IF  (INPT.GT.8)  GOTO  24 
20  CONTINUE 
C 

C  Input  of  the  structural  information.  N=#  of  modes,  NRB=#  of  rigid 
C  body  modes  (set  at  two  for  now),  NS=#  of  structural  modes, 

C  NC=f  of  control  surface  modes.  OMEGA-modal  frequencies 
C  GM-genera I i zed  mass,GK-genera I i zed  stiffness. 

C 

WRITE(*, ’ (A) ’)  *  Is  the  structural  information  in  a  file  ?’ 
READ(*,111)  Q1 
IF  (Ql.EQ.’Y’)  THEN 

WRITE (* , ’ (A) ’)  ’  Please  input  the  name  of  the  file.’ 

READ(* ,111)  STRUCT 

OPEN (1 , FILE=STRUCT , STATUS= ’ OLD  * ) 

READ (1,116)  N , NRB , NS , NC 
N2=2*N 

READ (1,117)  (GM(I) , 1=1 ,N) 

READ (1,117)  (OMEGA (I), 1=1, N) 

ELSE 

WRITE(*, ’ (A) ’)  ’  Do  you  wish  to  save  the  data  on  a  file  for’ 
WRITE(* ,  ’  (A)  ’)  ’  future  use?’ 

READ (*,m)  Q2 
IF  (Q2.EQ. ’Y’)  THEN 

WRITE (*, ’ (A) ’)  ’  What  is  the  name  of  the  new  file?’ 

READ (*,111)  STRUCT 

OPEN (1 , FILE=STRUCT , STATUS= ’ NEW ’ ) 

ENDIF 

WRITE(*, ’ (A) ’)  *  Input  total  f  of  modes  used  (max  20).’ 

READ(*,*)  N 
N2=2*N 

WRITE(* , ’ (A) ’)  *  Input  f  of  rigid  modes  used  (2).' 

READ(*,*)  NRB 

WRITE(* , * (A) ’)  ’  Input  |  of  structural  modes  used  (max  15).’ 
READ(*, *)  NS 

WRITE(* , ’ (A) ’)  ’  Input  f  of  control  surface  modes  used  (max  5).’ 
READ(*, *)  NC 

WRITE (* , ’ (A) ’)  ’  Input  the  generalized  masses  in  the  order  of’ 
WRITE(* , ’ (A) ’)  ’  rigid  modes  (plunge  then  pitch),  structural’ 
WRITE (* , ’ (A) ')  ’  modes,  and  control  surface  modes.’ 

READ(* , *)  (GM(I) ,I=1,N) 

WRITE(*, ’ (A) ’)  ’  Input  the  modal  frequencies  (in  Hz) in  the  order’ 
WRITE(* , ’ (A) ’)  ’  of  rigid  modes  (plunge  then  pitch),  structural  ’ 
WRITE (* , ’ (A) ’)  ’  modes,  and  control  surface  modes.’ 
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READ(*,*)  (OMEGA (I), 1=1, N) 

IF  (Q2.EQ. ’Y’)  THEN 
WRITE(1 ,116)  N , NRB , NS , NC 
WRITE(1 ,117)  (GM(I) ,I=1,N) 

WRITE(1 , 117)  (OMEGA (I), 1=1, N) 

END  IF 
ENDIF 
GOTO  24 

21  CONTINUE 
C 

C  Input  of  the  generalized  aerodynamic  forces.  It  is  in  the  format 
C  of  the  output  of  FASTOP  as  modified  by  Max  Blair,  where  the  two 
C  reduced  frequencies  are  combined  into  one  file. 

C  K1  is  the  first  reduced  frequency,  K2  is  the  second  reduced  frequency. 

C  GFl  are  the  generalized  forces  related  to  the  first  reduced  frequency, 

C  GF2  are  the  generalized  forces  related  to  K2.  K1  is  usually  0.0  and 

C  K2  is  a  small  reduced  frequency  (0.001). 

C 

WRITE (*, ’ (A) ’)  ’  What  is  the  name  of  the  file  for  the  generalized’ 
WRITE(*, ’ (A) ’)  ’  forces  (aerodynamics)?’ 

READ (*,111)  AERO 

OPEN (2 , FILE=AER0 , STATUS= ’ OLD ’ ) 

READ(2,111)  TITLE 
READ (2, 128)  N , XMACH 
READ (2, 127)  K1 

READ (2, 126)  (((F(I, J) ,G(I, J)) ,I=1,N) , J=1,N) 

DO  17  1=1, N 
DO  17  J=1,N 

GFl (I , J) =CMPLX (F (I , J) , G (I , J) ) 

17  CONTINUE 

READ (2, 127)  K2 

READ (2, 126)  ( ( (F (I , J) , G (I , J) ) , 1=1 , N) , J=1 , N) 

DO  18  1=1, N 
DO  18  J=1,N 

GF2 (I , J) =CMPLX (F (I , J) , G (I , J) ) 

18  CONTINUE 
8  CONTINUE 

GOTO  24 

22  CONTINUE 

Input  of  initial  conditions,  velocity,  density,  reference  area, 
reference  chord  (should  be  the  actual  A/C  area  £  chord). 

Also  the  appropriate  magnitudes  of  the  mode  shapes  are  input. 

WRITE (* , ’ (A) ’)  ’  Are  the  initial  conditions  in  a  file  ?’ 
READ(*,111)  Q3 
IF  (Q3.EQ. ’Y’)  THEN 

WRITE(* , ’ (A) ’)  ’  Please  input  the  name  of  the  file.’ 

READ(*, 111)  INCOND 

OPEN (3 , FILE=INCOND , STATUS= ’ OLD ’ ) 

READ(3, 118)  RHOO , AO , SREF , CBAR , BR 
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IF  (NRB.NE.O)  READ (3, 119)  ALPHAO , HO 
IF  (NS.NE.O)  READ (3 ,119)  (ZETAO(I) ,1=1, NS) 

IF  (NC.NE.O)  READ (3, 119)  (DELTAO(I) , 1=1 ,NC) 

ELSE 

WRITE(* , ’ (A) ’)  ’  Do  you  wish  to  save  the  initial  conditions  ’ 
WRITE(*, ’ (A) ’)  ’  on  a  file  for  future  use?’ 

READ(*,111)  Q4 
IF  (Q4.EQ. ’Y’)  THEN 

WRITE (* , ’ (A) ’)  ’  What  is  the  name  of  the  new  file?’ 

READ (*,111)  INCOND 

OPEN (3 , FILE=INCOND , STATUS=  * NEW ’ ) 

ENDIF 

WRITE (*, ’ (A) ’)  *  Please  input  the  following  initial  conditions. 
WRITE(* , ’ (A) ’)  ’  Density  (s I ugs/f t**3) : * 

READ(*,*)  RHOO 

WRITE(*, ’ (A) ’)  ’  Speed  of  sound  (ft/sec)  at  density  altitude:’ 
READ(*,*)  AO 

WRITE(*, ’ (A) ’)  ’  Reference  Area  (in  ft**2) : ’ 

READ(*, *)  SREF 

WRITE (*, ’ (A) ’)  ’  Reference  Chord  (MAC)  (in  ft):’ 

READ(*,*)  CBAR 

WRITE(*, ’ (A) ’)  ’  What  is  the  reference  semichord  used  in  the’ 
WRITE(*, ’ (A) ’)  ’  aerodynamics  routine  (from  k=bw/V)  in  ft.’ 

READ(*,*)  BR 
IF  (NRB.Eq.O)  GOTO  19 

WRITE(*, ’ (A) ’)  ’  Pitch  mode  shape  magnitude:’ 

READ(*, *)  ALPHAO 

WRITE (* , ’ (A) ’)  ’  Plunge  mode  shape  magnitude:’ 

READ(*,*)  HO 
19  CONTINUE 

IF  (NS.EQ.O)  GOTO  5 
DO  5  1=1, NS 
WRITE (*,108) I 
READ(*, *)  ZETAO(I) 

5  CONTINUE 

IF  (NC.EQ.O)  GOTO  6 
DO  6  1=1, NC 
WRITE (*,109) I 
READ(*,*)  DELTAO(I) 

6  CONTINUE 
IF(Q4.Eq. ’Y’)  THEN 

WRITE(3, 118)  RHOO, AO, SREF, CBAR, BR 
IF  (NRB.NE.O)  WRITE (3, 119)  ALPHAO, HO 
IF  (NS.NE.O)  WRITE (3, 119)  (ZETAO(I) ,1=1, NS) 

IF  (NC.NE.O)  WRITE (3, 119)  (DELTAO (I) , 1=1 , NC) 

ENDIF 
ENDIF 
GOTO  24 
26  CONTINUE 
C 

C  Manual  input  of  velocity  root  locus  data 


c 


WRITE(* , ’ (A) ’)  *  Do  you  wish  to  perform  a  velocity  root  locus?’ 
READ (*,100)  Q5 
IF(Q5.EQ. ’Y’)  THEN 

^  WRITE (* , ’ (A) ’)  ’  Are  the  values  on  file?’ 

READ (*,100)  Q6 
IF  (Q6 . EQ. ’Y ’)  THEN 

WRITE (*, ’ (A) ’)  ’  What  is  the  name  of  the  file?’ 

READ (*,100)  RTLCSF 

OPEN (7 , FILE=RTLCSF , STATUS= ’ OLD ’ ) 

\  READ (7, 121)  NVEL 

DO  16  1=1, NVEL 
READ (7, 122)  V(I) 

16  CONTINUE 
|  ELSE 


WRITE(* , ’ (A) ’)  ’  Do  you  want  to  save  this  data  for  future  use?’ 
READ(*,111)  Q7 
IF(Q7.EQ.’Y’)  THEN 

WRITE (* , ’ (A) ’)  ’  What  is  the  name  cf  the  new  file?’ 

READ (*,111)  RTLCSF 

OPEN (7 , FILE=RTLCSF , STATUS= ’ NEW ’ ) 

ENDIF 

WRITE (* , ’ (A) ’)  ’  How  many  velocities  tc  use  (max  20)?’ 

READ(*, *)  NVEL 

IF  (Q7.EQ. ’Y’)  WRITE (7, 101)  NVEL 

WRITE (* , ’ (A) ’)  ’  Please  input  the  velocities  (in  ft/sec).’ 
READ(* , *)  (V (I), 1=1, NVEL) 

IF  (Q7 . EQ. * Y * )  WRITE (7, 102)  (V (I) ,1=1, NVEL) 

ENDIF 
ENDIF 
GOTO  24 
25  CONTINUE 
C 

C  Manual  input  of  the  forcing  function  matrix 
C 

WRITE(* , ’ (A) ’)  ’  Is  the  forcing  function  matrix  on  file?’ 
READ(*,111)  Q8 
IF(Q8.EQ. ’Y’)  THEN 

WRITE (* , * (A) ’)  ’  Please  input  the  name  of  the  file.’ 

READ(* , 111)  INPTM 

OPEN (11 , FILE=INPTM , STATUS= ’ OLD ’ ) 

READ(11 , 121)  NINP 

READ (11,123)  ( (QMAT (I , J) , J=1 , NINP) , 1=1 , N) 

ELSE 

WRITE(* ,  ’  (A) ’)  ’  Do  you  wish  to  save  the  data  to  file?’ 

READ(* , 111)  Q9 
IF  (Q9.EQ. ’ Y * )  THEN 

WRITE (* , ' (A) ’)  ’  What  is  the  name  of  the  new  file?’ 

READ(*, 111)  INPTM 

OPEN (1 1 , FILE=INPTM , STATUS= ’ NEW ’ ) 

ENDIF 
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WRITE (*, * (A) ’)  *  How  many  inputs  are  there?’ 

READ(*,*)  NINP 
WRITE (*,124)  N,NINP 
00  30  1=1,  N 
WRITE (*,125)  I 

READ(*,*)  (QMAT (I , J) , J=1 , NINP) 

30  CONTINUE 

IF  (Q9.EQ. *Y’)  THEN 
WRITE(11,121)  NINP 

WRITE (11 ,123)  ( (QMAT (I , J) , J=1 , NINP) , 1=1 , N) 

ENDIF 
ENDIF 
GOTO  24 
27  CONTINUE 

Manual  input  of  the  C  matrix 

WRITE (*, ’ (A) ’)  ’  Is  the  C  matrix  on  file?’ 

READ(*,111)  Q10 
IF(Q10.EQ. ’Y’)  THEN 

WRITE(*, ’ (A) ’)  ’  Please  input  the  name  of  the  file.’ 

READ(*,111)  OUTPM 

OPEN (12 , FILE=OUTPM , STATUS* *  OLD ’ ) 

READ (12, 121)  NOUT 

READ (12 , 123)  ( (CMAT (I , J) , J=1 , N2) , 1=1 , NOUT) 

ELSE 

WRITE(*, ’ (A) ’)  ’  Do  you  wish  to  save  the  data  to  file?’ 
READ(*,111)  Qll 
IF  (Qll.EQ.’Y’)  THEN 

WRITE(*, ’ (A) ’)  ’  What  is  the  name  of  the  new  file?’ 

READ (*,111)  OUTPM 

OPEN (12 , FILE=OUTPM , STATUS* ’ NEW ’ ) 

ENDIF 

WRITE(*, ’ (A) ’)  ’  How  many  outputs  are  there?’ 

READ(*,*)  NOUT 
WRITE (*,129)  N0UT,N2 
DO  34  1=1, NOUT 
WRITE (*,125)  I 
READ(*,*)  (CMAT (I , J) , J=1 , N2) 

34  CONTINUE 

IF  (Qll.EQ.’Y’)  THEN 
WRITE (12, 121)  NOUT 

WRITE (12 , 123)  ( (CMAT (I , J) , J=1 , N2) , 1=1 , NOUT) 

ENDIF 
ENDIF 
GOTO  24 
23  CONTINUE 
DO  9  1=1, N 

GK (I) =GM (I) * ( (OMEGA (I) *2 . *PI) **2) 

9  CONTINUE 
VO=XMACH*AO 


WRITE (4, 100)  N 
WRITE (4, 101)  NRB 
WRITE (4, 102)  NS 
WRITE (4, 103)  NC 
WRITE (4, 104) 

DO  1  1=1, N 

WRITE (4, 105)  I , I , GM (I) 

1  CONTINUE 
WRITE (4, 106) 

DO  2  1=1, N 

WRITE (4, 105)  I , I , GK (I) 

2  CONTINUE 
C 

C  Call  to  the  routine  that  non-dimensional izes  FASTOP  aero  forces 
C 

CALL  FASTCHG 
IF  (NAERO.EQ.O)  GOTO  4 
WRITE (4, 114)  XMACH 
WRITE (4, 107)  K1 
DO  3  1=1, N 
WRITE (4, 125)  I 

WRITE (4 ,112)  ( (REAL (GFl (I , J) ) , AIMAG (GFl (I , J) ) ) , J=1 , N) 

3  CONTINUE 
WRITE (4, 107)  K2 
DO  4  1=1, N 

WRITE (4, 125)  I 

WRITE (4,112)  ( (REAL (GF2 (I , J) ) , AIMAG (GF2 (I , J) ) ) , J=1 , N) 

4  CONTINUE 

WRITE(4, ’ (/A) ’)  ’  The  initial  conditions  are:’ 

WRITE (4, 115)  RHOO , SREF , CBAR , BR 
WRITE (4, ’  (A)  ’)  *1’ 

100  FORMAT (/,’  The  total  number  of  modes  are:’, 13) 

101  FORMAT(’  The  total  number  of  rigid  body  modes  are: ’,13) 

102  FORMAT(’  The  total  number  of  elastic  modes  are: ’,13) 

103  FORMAT (’  The  total  number  of  control  surface  modes  are: ’,13) 

104  FORMAT (/,’  The  generalized  mass  matrix  is’,/,’  RowCol  Gen  Mass’) 

105  FORMAT (IX, 213, E12. 4) 

106  FORMAT (/,’  The  generalized  stiffness  matrix  is’, 

0/,’  RowCol  Gen  Stiff’) 

107  FORMAT (/,’  The  generalized  force  matrix  for  reduced  frequency  of’, 
0F7 . 4 , ’  is’) 

108  FORMAT (’  Structural  mode ’,12,’  magnitude:’) 

109  FORMAT (’  Control  surface  mode ’,12,’  magnitude:’) 

110  FORMAT (IX, 14 , IX, 14 , 2E15 . 7) 

111  FORMAT (A) 

112  FORMAT (3 (IX , 2E12 . 4 , ’ i ’ ) ) 

113  FORMAT (15X, El 5. 7) 

114  FORMAT(/,’  For  a  Mach  number  of’,F6.3) 

115  FORMAT(’  Densi ty=’ ,Fll .5, ’SI ugs/ft**3,  Reference  Areas’, 

0E11.4,’  ft**2’ ,/, 

0’  Chords’ , Ell .4, *  feet,  and  an  aerodynamic  reference’, 


O’  sem i -chords’ , Ell .4, *  feet.’) 

116  FORMAT (415) 

117  FORMAT (4E15. 7) 

118  FORMAT (F10.7,F10.2,3F10.5) 

119  FORMAT (6F10. 5) 

120  FORMAT (’  MAC  ready  for  input',//, 

0’  Enter  (#) :  ’,//, 

O'  1.  Automatic  input  of  items  (2) -(7)’,/ 

0’  2.  Read  in  structural  data’,/ 

0’  3.  Read  in  aerodynamic  data’,/ 

0’  4.  Read  in  initial  conditions’,/ 

0’  5.  Read  in  root  locus  values’,/ 

O'  6.  Read  in  forcing  function  matrix',/ 

0’  7.  Read  in  the  C  matrix’,/ 

0’  8.  Completion  of  input  data’,//) 

121  FORMAT (15) 

122  FORMAT (F10. 3) 

123  FORMAT (6E15. 7) 

124  FORMAT (  ’  For  n  modes,  there  are  n  modes  with  m  controls’,/, 
0’  The  modes  will  be  set  up  as  follows:’,/, 

0’  Rigid  modes’ ,/, 

0’  Structural  modes’,/, 

0’  Control  surface  modes’,/, 

0’  The  matrix  will  have ’,13,’  rows  and  ’,13,'  columns’,/, 

0’  Please  input  the  matrix  by  row.’//) 

125  FORMAT (’  Row ’,13) 

126  FORMAT (8E15. 7) 

127  FORMAT (E15. 7) 

128  FORMAT (15 , 5X , E15 . 7) 

129  F0RMAT(’  The  matrix  will  have  ’,13,’  rows  and', 13,'  columns’, 
0’  Please  input  the  matrix  by  row.’//) 

RETURN 

END 


ooooonooooooooooooooonooooooooooooooooooooooooo 


SUBROUTINE  MATSAV  (  LUNIT,  NAME,  NR,  M,  N,  IMG, 
$  XREAL,  XIMAG,  FORMT  ) 


(TM) 

MATRIX  V5.0 
X 

(C)  COPYRIGHT  INTEGRATED  SYSTEMS,  INC.  1985 
PALO  ALTO,  CALIFORNIA 

All  Rights  Reserved 

This  software  may  not  be  copied  or  altered 
without  the  express  written  consent  of 
Integrated  Systems,  Inc. 


MATSAV  writes  a  matrix  to  a  file  in  a  format  suitable  for  the 
MATRIXx  LOAD  operation. 


Pa ram. 

1  Type 

|  On  input- 

|  On  output- 

LUNIT 

1 

|  INTEGER 

i 

|  Fortran  logical  unit  number. 

1 

|  unchanged. 

1 

NAME 

|  CHARACTER* (*) 
j  (maximum 
j  1 ength  10) 

I 

|  Name  of  the  matrix.  One  al- 
j  phabetic  followed  by  up  to  9 
j  alphanumeric  characters. 

1 

1 

|  unchanged. 

1 

1 

1 

NR 

|  INTEGER 

1 

1 

1 

1 

1 

|  Row-dimension  in  the 
j  defining  dimension  or  type 
j  statement  i n  the  calling 
j  program.  NR  must  be  greater 
j  than  or  equal  to  M. 

i 

|  unchanged. 

1 

1 

1 

1 

1 

M 

|  INTEGER 

|  Number  of  rows  of  the  matrix 

1 

|  unchanged. 

1 

N 

|  INTEGER 

1 

|  Number  of  columns  of  the 
j  matrix. 

|  unchanged. 

1 

1 

IMG 

|  INTEGER 

1 

1 

i 

1 

|  If  IMG  =  0,  the  imaginary 
j  part  (XIMAG)  is  assumed  to 
j  be  zero  and  is  not  saved. 

i 

1 

|  unchanged. 

1 

1 

l 

XREAL 

|  DOUBLE 

1 

|  Real  part  of  the  matrix  to 

1 

|  unchanged. 

or>or>r>oor>or»oooor»or>r>or>oooor>ooor»r»r>or»r>r>r» 


XIMAG 


FORMT 


PRECISION 


DOUBLE 

PRECISION 


CHARACTER* (*) 
(maximum 
length  20) 


be  saved. 


Imaginary  part  of  the  matrix 
to  be  saved. 


String  containing  the  For¬ 
tran  format  to  be  used  for 
writing  the  elements  of  the 
matr i x . 


unchanged . 


unchanged. 


Example:  The  following  Fortran  program  generates  an  elementary 
matrix  in  X  and  writes  it  to  Fortran  unit  1.  Assume 
that  unit  1  has  been  preallocated  as  file  (data  set) 
TEST. 


DOUBLE  PRECISION  X(20,3),  DUMMY 
DO  200  J=l,3 
DO  100  1=1,10 
X(I, J)=0.0D0 
100  CONTINUE 
X(J, J)=l .000 
200  CONTINUE 

CALL  MATSAV (  1,  ’AMATRIX’,  20,  10,  3,  0, 
S  X,  DUMMY,  * (1P2E24.15) ’  ) 

STOP 
END 


After  this  program  runs,  invoke  MATRIXx  and  type: 


<>  LOAD  ’TEST’ 


This  will  put  X  on  the  stack  as  stack-variab I  e-name  AMATRIX. 


INTEGER  LUNIT,  M,  N,  NR,  IMG 

CHARACTER* (*)  NAME,  FORMT 

DOUBLE  PRECISION  XREAL(NR,1) ,  XIMAG(NR,1) 

CHARACTER  NAM* 10,  FORM* 20 


|  write  header  record. 


NAM=NAME 

FORM=FQRMT 

WRITE (LUNIT, ’(A10, 315, A20) ’)  NAM,M,N, IMG, FORM 


V.V.V.  v-Avv.y 


•V.< 


*V<' 


l.v 


ks!»s 


c 


|  write  real -part  of  the  matrix 


WRITE (LUNIT, FORM)  ( (XREAL (I , J) , 1=1 , M) , J=1 , N) 


SUBROUTINE  MINV 
C 

C  This  subroutine  inverts  the  X  matrix  using  the  IMSL  routine 
C  LINV1F. 

C 

REAL  X,XINV,WKAREA 

DIMENSION  X(20,20)  ,XIM(20,20)  ,WKAREA(20) 

COMMON/INV/X , XINV 
INTEGER  N,NRB,NS,NC,IER 
COMMON/INDEX/N , NRB, NS , NC 
INTEGER  NAERO , NMAT , NINV , NAMAT , NEIG 
COMMON/PRT /NAERO , NMAT, NINV , NAMAT, NEIG 
C 

C  Call  to  IMSL  matrix  inversion  routine 
C 

CALL  LINV1F (X , N , 20 , XINV , 4 , WKAREA , IER) 

IF  (IER.EQ.34)  THEN 
WRITE (*,403) 

GOTO  5 
ENDIF 

IF  (IER.NE.O)  THEN 
WRITE (*,400)  IER 
STOP 
ENDIF 

5  CONTINUE 
IF(NINV.EQ.O)  GOTO  2 
WRITE (4, * (/A) ’)  ’  X  inverse  is’ 

DO  2  1*1, N 
WRITE (4, 401)  I 
DO  3  J=1,N,6 

WRITE (4, 402)  (XINV (I , K) , K=J , J*5) 

3  CONTINUE 
2  CONTINUE 

400  FORMAT(///,’  ****  There  was  an  error  in  inverting  X,  f  ’,13, 

0’  ****’,///) 

401  FORMAT (’  ROW ’,13) 

402  FORMAT (6E12. 4) 

403  FORMAT(///,’  The  inversion  did  not  meet  error  criterion’,///) 
1  CONTINUE 

RETURN 

END 


o  o  nnnn  r>  o  o  o 


SUBROUTINE  MMULT 


This  subroutine  multiplies  Xinv*Y  and  Xinv*Z  using  the 
IMSL  routine  VMULFF 

REAL  Y ,  Z , XINV , XINVY , XINVZ 

DIMENSION  X (20 , 20) , Y (20 , 20) , Z (20 , 20) , XINV (20 , 20) , 

OXINVY (20 , 20) , XINVZ (20 , 20) 

COMMON/MMUL/Y , Z , XINVY , XINVZ 

COMMON/INV/X,XINV 

INTEGER  N,NRB,NS,NC,IER 

COMMON/INDEX/N , NRB , NS , NC 

INTEGER  NAERO , NMAT , NINV , NAMAT , NEIG 

COMMON/PRT/NAERO , NMAT , NINV , NAMAT , NEIG 

Call  to  IMSL  matrix  multiplication  routine  to  multiply 
XINV*Y 

CALL  VMULFFCXINV.Y.N.N.N^O^O.XINVY^O.IER) 

IF  (IER.NE.O)  THEN 
WRITE (*,500)  IER 
STOP 
ENDIF 

IF (NINV . EQ . 0)  GOTO  2 
WRITE (4, ’(/A)’)  »  XINV*Y  is’ 

DO  2  1=1, N 
WRITE (4, 501)  I 
DO  3  J=1,N,6 

WRITE (4, 502)  (XINVY (I, K) ,K=J, J*5) 

3  CONTINUE 
2  CONTINUE 

IER=0 

Call  to  IMSL  matrix  multiplication  routine  to  multiply 
XINV*Z 

CALL  VMULFF (XINV, Z,N,N,N, 20, 20, XINVZ, 20, IER) 

IF  (IER.NE.O)  THEN 
WRITE (*,500)  IER 
GOTO  1 
ENDIF 

IF(NINV.EQ.O)  GOTO  4 
WRITE (4, ’ (/A) ’)  ’  XINV*Z  is’ 

DO  4  1=1, N 
WRITE (4, 501)  I 
DO  5  J=1 ,N,6 

WRITE (4, 502)  (XINVZ(I,K) ,K=J, J*5) 

5  CONTINUE 

4  CONTINUE 

500  FORMAT(///,’  ****There  was  an  error  in  the  matrix  multiply,  f  ’ 


0,13/  **.-////) 

501  FORMAT (’  ^JW’,13) 

502  FORMAT  (6E:  '1 . 4) 

1  CONTINUE 

RETURN 

END 


nnn  o  oo  r>  r>  r> 


SUBROUTINE  RLPLOT (FL , FH , RL , RH , TITLE) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*4  XPOS , YPOS 
CHARACTER *4  NUM 
CHARACTER*23  DT 
CHARACTER *40  TITLE 
LOGICAL  PLOT 

REAL  VAL(13)/.01, .1, .2, .5,1. ,2. ,5. ,10. ,20. ,50. ,100. ,200. ,500./ 

C 

REWIND  90  i 

ZER0=O . DO 
0NE=1 . DO 
WRITE (*,1010) 

READ(*,*) 

CALL  JBEGIN 

C  CLEAR  SCREEN  WITH  PLOTlO  COMMANDS 
CALL  INITT (960) 

CALL  FINITT (0,2800) 

SETUP  GRAPHICS 

PLOT=. FALSE. 

'  IF (.NOT.  PLOT)  THEN 
CALL  JDINIT (1) 

CALL  JDEVON(l) 

ELSE 

CALL  JDINIT (2) 

CALL  JDEVON (2) 

END  IF 

CALL  JWINDO (-4 . , 12 . , -3 . , 15 . ) 

CALL  JOPEN 
CALL  JSIZE(.27, .35) 

IF (PLOT)  CALL  JF0NT(5) 

IF(. NOT. PLOT)  CALL  JFONT(l) 

DRAW  BOX  AND  SET  UP  AXES 

CALL  JMOVEro.,0.) 

CALL  JDRAW(10. ,0.) 

CALL  JDRAW(10. ,10.) 

CALL  JDRAW(0. ,10.) 

CALL  JDRAW (0 . , 0 . ) 

XRANGE=RH-RL 
XRAT=XRANGE/10. 

YRANGE=FH-FL 
YRAT=YRANGE/10. 

C 

DO  50  1=1,13 


IF (XRANGE . GE . VAL (I) )  XSTEP=VAL (I)/5. 

50  IF (YRANGE . GE . VAL (I) )  YSTEP=VAL(I)/5. 

C 

XP0S=SNGL (-RL/XRAT) 

CALL  JMOVE (XPOS , 0 . ) 

CALL  JDRAW(XPOS,10.) 

C 

C  DRAW  TIC  MARKS  AND  NUMBERS 
C 

X=RL 

60  XPOS=SNGL (X/XRAT -RL/XRAT) 

IF  (XPOS. GT. 10.) GOTO  65 
CALL  JMOVE (XPOS ,0.) 

CALL  JRDRAW (0 . ,  - .  3) 

C 

ENCODE (4, 1020, NUM)X 
CALL  JJUST (2,3) 

CALL  JMOVE (XPOS, -.5) 

CALL  JHSTRG (NUM) 

C 

X=X+XSTEP 
GOTO  60 
C 

65  Y=FL 

70  YPOS=SNGL (Y/YRAT -FL/YRAT) 

IF(YPOS.GT.10.)  GOTO  75 
ENCODE (4, 1020, NUM) Y 
CALL  JJUST (3,2) 

CALL  JMOVE (0. ,YPOS) 

CALL  JRDRAW (-.3,0.) 

CALL  JM0VE(-.5,YP0S) 

CALL  JHSTRG (NUM) 

Y=Y+YSTEP 
GOTO  70 
C 

75  CALL  JSIZE(.4, .45) 

CALL  JJUST (2, 3) 

CALL  JBASE(1 . ,0. ,0.) 

CALL  JMOVE (5., -2.) 

CALL  JHTEXT (42 , 42HR [BLC] EAL  [ELC] C [BLC] OMPONENT  (1/SEC) [ELC]) 
CALL  JBASE (0 . , 1 . , 0 . ) 

CALL  JJUST (2,1) 

CALL  JMOVE (-2., 5.) 

CALL  JHTEXT (44, 44HI [BLC] MAGINARY  [ELC] C [BLC] OMPONENT  (HZ) [ELC]) 
C 

C  PRINT  TITLE  AND  DATE  AT  TOP  OF  PLOT 
C 

C  CALL  JBASE (1 . ,0. ,0.) 

C  CALL  JJUST (1,1) 

C  CALL  JMOVE (-3., 12.) 


CALL  JHSTRG (TITLE) 

C  CALL  LIB$DATE_TIME (DT) 

C  CALL  JM0VE(-3. ,11.) 

C  CALL  JHSTRG (DT) 

C 

80  READ (90 , END=999)  XH,YH,RF 

IF (XH .GT.RH.OR.XH.LT. RL)  GOTO  80 
IF (YH .GT.FH.OR.YH.LT. FL)  GOTO  80 
XPOS=SNGL (XH/XRAT -RL/XRAT) 

YPOS=SNGL (YH/YRAT-FL/YRAT) 

CALL  JMOVE (XPOS , YPOS) 

C 

C  DRAW  A  LITTLE  SQUARE 
C 

CALL  JRDRAW (.05, .05) 

CALL  JRDRAW(- .1,0.) 

CALL  JRDRAW(0. ,  - .  1) 

CALL  JRDRAW (.1,0.) 

CALL  JRDRAW (0. , .1) 

GOTO  80 
C 

C  SHUT  EVERYTHING  OFF 
C 

999  CALL  JCLOSE 

IF (.NOT. PLOT)  THEN 
CALL  JPAUSE(l) 

CALL  JDEVOF(l) 

CALL  JDEND(l) 

ELSE 

CALL  JPAUSE(2) 

CALL  JDEVOF (2) 

CALL  JDEND(2) 

END  IF 
C 

C  SEND  QMS  FILE  TO  THE  PRINTER 
C 

C  IF (PLOT) THEN 

C  WRITE (*,1040) 

C  ISPAWN=LIBSSPAWN ( ’ Q [MAX . ASE . ADAM . DI3SPAWN] MAXRL . COM ’ ) 

C  IF (.NOT. I SPAWN)  CALL  LIBSSIGNAL(XVAL(ISPAWN)) 

C  PLOT=. FALSE. 

C  END  IF 

C 

C  MENU 

C  CLEAR  SCREEN  WITH  PLOTIO  COMMANDS 
210  CALL  INITT (960) 

CALL  FINITT (0,2800) 

WRITE (*,1030) 

READ(*,*)IGOTO 

IF (IGOTO . LT . 1 . OR . IGOTO . GT . 3)  GOTO  210 
GOTO  (300, 400, 500), IGOTO 


C  RETURN  TO  ADAM 
300  CALL  JEND 

WRITE (*,1050) 

READ(*,*) 

RETURN 

C 

C  CHANGE  PLOT  LIMITS 
C 

400  WRITE (*,1060)  RL,RH,FL,FH 
WRITE (*,1070) 

READ (*, 1065) RL 
WRITE (*,1080) 

READ (*, 1065) RH 
WRITE (*,1090) 

READ (*, 1065) FL 
WRITE (*,1100) 

READ (*, 1065) FH 
REWIND  90 
GOTO  30 
C 

C  MAKE  HARDCOPY  OF  CURRENT  PLOT 
500  PLOT=.TRUE. 

REWIND  90 
GOTO  30 


C 

1010 

1020 

1030 


FORMAT (’  SWITCH  TO  GRAPHICS  MODE  AND  HIT  RETURN’) 

FORMAT (F4.0) 

FORMAT (///’  PLOT  ROUTINE  MENU’//’  1-RETURN  TO  MAC’//’  2-CHANGE 
PLOT  LIMITS’//’  3-MAKE  A  HARDCOPY  OF  CURRENT  PLOT’//) 

FORMAT (’  SUBPROCES  SPAWNED.  FILE  ON  ITS  WAY  TO  PLOTTER’) 

FORMAT (’  RETURN  TERMINAL  TO  STANDARD  MODE  AND  HIT  RETURN  KEY’) 
FORMAT ( ’  XMIN=’,F8.4/’  XMAX= ’ , F8 . 4/ ’  YMIN=’,F8.4/’  YMAX=’,F8.4) 
FORMAT (F9. 5) 

FORMAT ( ’  NEW  XMIN  ?’/) 

FORMAT (’  NEW  XMAX  ?’/) 

FORMAT (’  NEW  YMIN  ?’/) 

FORMAT (’  NEW  YMAX  ?’/) 

FORMAT (’  NEW  ROW  ?’/) 

FORMAT (’  NEW  COL  ?’/) 

FORMAT (13) 

END 
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SUBROUTINE  STABDER (NVEL.V) 


This  subroutine  computes  the  dimensional,  and  some  non-dimensional 
stability  derivatives  for  use  in  developing  the  reduced  order  model. 

CHARACTER *1  QCHANGE 
INTEGER  N,NRB,NS,NC 
COMMON/INDEX/N , NRB , NS , NC 
REAL  K1,K2 
COMPLEX  GF1.GF2 

DIMENSION  GF1 (20 , 20) , GF2 (20 , 20) 

C0MM0N/AER0/K1 , K2 , GFl , GF2 , PI 

REAL  VO , RHOO , SREF , CBAR , ALPHAO , HO , ZETAO , DELTAO , BR 

DIMENSION  ZETAO (15), DELTAO (5) 

COMMON/IC/VO , RHOO , SREF , CBAR , ALPHAO , HO , ZETAO , DELTAO , BR 
REAL  ZALPHA , MALPHA , ZALPHDT , MALPHDT ,ZQ,MQ, 

OFZALPHA , FZALPDTN , FZETQ , ZZETA , MZETA , 

OZZETADT ,  MZETADT ,  ZZE7DDT ,  MZETDDT , 

OFZETZET , FZETZDT , FZEZDDT , 

OFZETDEL , FZETDDT , FZEDDDT , 

OFDELZET, FDELZDT, FDEZDDT, 

OZDELTA , MDELTA , FD ALPHA , FDALPDT , FDETQ , 

OZDELTDT , MDELTDT , ZDELDDT , MDELDDT , 

OFDELDEL , FDELDDT , FDEDDDT , VEL 

DIMENSION  FZALPHA(15) , FZALPDT (15) , FZETQ (15) , ZZETA (15) , MZETA (15) , 
OZZETADT (15) , MZETADT (15) ,ZZETDDT(15) , MZETDDT (15) , 

OFZETZET (15, 15) , FZETZDT (15, 15) .FZEZDDT (15, 15) , 

OFZETDEL (15, 5) .FZETDDT (15, 5) .FZEDDDT (15, 5) , 

OFDELZET (5,15), FDELZDT (5 , 15) .FDEZDDT (5, 15) , 

OZDELTA (5) , MDELTA (5) , FDALPHA (5) , FDALPDT (5) , FDETQ (5) , 

OZDELTDT (5) .MDELTDT (5) .ZDELDDT (5) .MDELDDT (5) , 

OFDELDEL (5,5), FDELDDT (5,5), FDEDDDT (5 , 5) 

COMMON/DERIV/ZALPHA , MALPHA , ZALPHDT , MALPHDT , ZQ , MQ , FZALPHA , FZALPDT , 
OFZETQ .ZZETA , MZETA , ZZETADT , MZETADT , ZZETDDT , MZETDDT , FZETZET , 
OFZETZDT, FZEZDDT, FZETDEL, FZETDDT, FZEDDDT, FDELZET, FDELZDT, FDEZDDT, 
OZDELTA , MDELTA , FDALPHA , FDALPDT , FDETQ , 

OZDELTDT , MDELTDT , ZDELDDT , MDELDDT , FDELDEL , FDELDDT , FDEDDDT , VEL 
REAL  CMALPHA , CZALPHA , CMALPDT , CZALPDT , CMQ , CZQ , CMDELTA , 
OCZDELTA.MODE 

DIMENSION  MODE (20) .CMDELTA (5) ,CZDELTA(5) 

If  a  velocity  root  locus  is  desired,  then  the  velocity  takes  on 
the  velocities  desired,  otherwise  the  initial  velocity  is  used. 

IF  (NVEL.EQ.O)  THEN 
VEL=VO 
ELSE 
va=v 

ENDIF 

WRITE (*,203)  VEL, RHOO 
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WRITE (4, 203)  VEL,RHOO 
QBAR=0 . 5*RHOO*VEL**2 
IF  (Va.EQ.O.O)  THEN 
QBARsO. 0 
VEL=0.001 
ENDIF 

IF  (NRB.EQ.O)  GOTO  29 
MODE(1)=HO 
MODE (2) sALPHAO 
29  CONTINUE 

IF  (NS.EQ.O)  GOTO  1 
DO  1  Isl,NS 
J=I*NRB 

MODE(J)=ZETAO(I) 

1  CONTINUE 

IF  (NC.EQ.O)  GOTO  2 
DO  2  1=1, NC 
JsI+NRB+NS 
MODE  (  J)  sDELTAO  (I) 

2  CONTINUE 

Computation  of  ail  the  dimensional  derivatives  is  accomplished 
below.  They  are  based  on  the  paper  by  Rodden,  AIAA 
and  the  formula  for  the  generalized  forces  from  FASTOP 

IF (NRB.EQ.O)  GOTO  20 

ZALPHA=REAL (GFl (1,2)) *QBAR*SREF/ (MODE (1) *MODE (2) ) 

MALPHA=REAL (GFl (2,2)) *QBAR*SREF*CBAR/ (MODE (2) **2) 
ZALPHDTs-REAL (GF2 (1,1))* (SREF*QBAR*BR**2/ (K2**2*VEL) ) / 
0(M0DE(1)**2) 

MALPHDTs-REAL (GF2 (2,1))* (SREF*CBAR*QBAR*BR**2/ (K2**2*VEL) ) / 
0 (MODE (1) *MODE (2) ) 

ZQ=AIMAG (GF2 (1,2))* (SREF*QBAR*BR/ (K2*VEL) ) / (MODE (1) *MODE (2) ) 
OZALPHDT 

MQ=AIMAG(GF2  (2,2))  *  (SREF*CBAR*QBAR*BR/  (K2*VEL)  )  /  (MODE  (2)  **2) 
OMALPHDT 
20  CONTINUE 

IF  (NS.EQ.O)  GOTO  22 
DO  3  1=1, NS 
JsI+NRB 

IF  (NRB.EQ.O)  GOTO  21 

FZALPHA  (I)  sREAL  (GFl  ( J ,  2) )  *QBAR*SREF/  (MODE  ( J)  *MODE  (2) ) 
FZALPDT(I)  =-REAL  (GF2  ( J ,  1) )  *  (SREF*QBAR*BR**2/  (K2**2*VEL)  )  / 

0  (MODE(l)  *MODE(J)  ) 

FZETQ (I) =AIMAG (GF2 ( J , 2) ) * (SREF*QBAR*BR/ (K2* VEL) ) / (MODE ( J) * 
QMODE (2) ) -FZALPDT (I) 

ZZETA (I) sREAL (GFl (1 , J) ) *QBAR*SREF/ (MODE (1 ) *MODE ( J) ) 

MZETA (I) =REAL (GFl (2 , J) ) *QBAR*SREF*CBAR/ (MODE (2) *MODE ( J) ) 
ZZETADT (I)=AIMAG (GF2 (1 , J) ) * (SREF*QBAR*BR/ (K2*VEL) ) / 

0 (MODE ( 1 ) *MODE ( J) ) 

MZETADT (I) sAIMAG (GF2 (2 , J) ) * (SREF*CBAR*QBAR*BR/ (K2*VEL) ) / 


0 (MODE (2) .MODE (J)) 

ZZETDDT(I)=- (REAL (GF2 (1 ,  J)). (SREF.QBAR/ (MODE (1) .MODE ( J) ) ) 
O-ZZETA (I) ) / ( (K2.VEL/BR) ..2) 

MZETDDT (I) =- (REAL (GF2 (2 ,  J) ) * (SREF.CBAR.QBAR/ 

0 (MODE (2) .MODE ( J) ) ) -MZETA (I) ) / ( (K2.VEL/BR) *.2) 

21  CONTINUE 

DO  4  K=1 , NS 

l=k+nrb 

FZETZET (I , K)=REAL (GFl ( J ,  L) ) .QBAR.SREF/ (MODE(J) .MODE (L) ) 

FZETZDT (I , K) =AIMAG (GF2 ( J , L) ) * (SREF.QBAR.BR/ (K2.VEL) ) / 

0 (MODE (J) .MODE (L)) 

FZEZDDT (I , K) =- (REAL (GF2 ( J , L) ) * (SREF.QBAR / (MODE ( J) .MODE (L) ) ) - 
OFZETZET (I ,K))/( (K2.VEL/BR) **2) 

4  CONTINUE 

IF  (NC.EQ.O)  GOTO  5 
DO  5  K=1 , NC 

l=k+nrb+ns 

FDELZET (K , I) =RFAL (GFl (L,  J) ) .QBAR.SREF.CBAR/ (MODE(L) *MODE(J) ) 

FZETDEL (I ,  K) =REAL (GFl (J , L) ) .QBAR.SREF/ (MODE(J) *MODE(L)) 

FDELZDT (K , I) =AIMAG (GF2 (L , J) ) * (SREF.CBAR.QBAR.BR/ (K2.VEL) ) / 

0 (MODE (L) .MODE (J)) 

FZETDDT (I , K) =AIMAG (GF2(J, L) ) . (SREF.QBAR.BR/ (K2*VEL) ) / 

0 (MODE (J) .MODE (L)) 

FDEZDDT (K , I) =- (REAL (GF2 (L, J) ) . (SREF.CBAR.QBAR/ (MODE (L) * 

OMODE ( J) ) ) -FDELZET (K , I) ) / ( (K2.VEL/BR) ..2) 

FZEDDDT (I , K) =- (REAL (GF2 (J,L))* (SREF.QBAR/ (MODE ( J) .MODE (L) ) ) - 
OFZETDEL (I ,  K)  )  /  ( (K2.VEL/BR) **2) 

5  CONTINUE 

3  CONTINUE 

22  CONTINUE 

IF  (NC.EQ.O)  GOTO  23 
DO  6  1=1, NC 
JsI.NRB.NS 

IF  (NRB.EQ.O)  GOTO  25 

ZDELTA (I) =REAL (GFl (1 , J)) .QBAR.SREF/ (MODE(l) .MODE(J) ) 
MDELTA(I)=REAL (GFl (2, J) ) *QBAR*SREF*CBAR/ (MODE (2) .MODE(J)) 
FDALPHA (I) =REAL (GFl ( J , 2) ) .QBAR.SREF.CBAR/ (MODE ( J) .MODE (2) ) 
FDALPDT (I) =-REAL (GF2 ( J , 1) ) . (SREF*CBAR*QBAR*BR**2/ 

0(K2*.2.VEL) ) / (MODE(l) *MODE(J) ) 

FDETQ(I) =AIMAG (GF2 ( J , 2) ) . (SREF.CBAR.QBAR.BR/ (K2*VEL) ) / 

0 (MODE ( J) .MODE (2) ) -FZALPDT (I) 

ZDELTDT (I) sAIMAG (GF2 ( 1 , J))* (SREF.QBAR.BR/ (K2.VEL) ) / 

0 (MODE (1). MODE (J)) 

MDB.TDT (I) =AIMAG (GF2 (2, J)) * (SREF.CBAR.QBAR.BR/ (K2.VEL))/ 

0 (MODE (2) .MODE(J)) 

ZDELDDT (I) =- (REAL (GF2 (1 , J) ) . (SREF.QBAR/ (MODE(l) .MODE ( J) ) ) - 
OZDELTA (I) ) / ( (K2.VEL/BR) ..2) 

MDELDDT(I) =- (REAL (GF2 (2 , J) ) . (SREF.CBAR.QBAR/ 

Q (MODE (2) .MODE ( J) ) ) -MDELTA(I) ) / ( (K2. VEL/BR) * .2) 

25  CONTINUE 

IF  (NC.EQ.O)  GOTO  24 


DO  6  K=1,NC 
L=K+2*NS 

FDELDEL (I , K) =REAL (GFl ( J , L) ) * QBAR* SREF *CBAR/ (MODE ( J) *MODE (L) ) 
FDELDDT (I , K) =AIMAG (GF2(J, L) ) * (SREF*CBAR*QBAR*BR/ (K2*VEL) ) / 

0 (MODE (J) .MODE (L)) 

FDEDDDT (I , K) =- (REAL (GF2 ( J , L) ) * (SREF.CBAR.QBAR/ 

0(M0DE(J) *MODE(L) ) ) -FDELDEL (I, K) ) / ( (K2.VEL/BR) **2) 

24  CONTINUE 

6  CONTINUE 
23  CONTINUE 

C 

C  Computation  of  the  non-dimensional  stability  derivatives  of  interest 
C  from  the  dimensional  derivatives.  These  can  be  used  for  comparison 
C  against  another  aerodynamic  code  to  check  for  accuracy. 

C 

IF  (NRB.EQ.O)  GOTO  26 

CMALPHA=MALPHA / (0 . 5*RHOO.VEL*.2*SREF*CBAR) 

CZALPHA=Z ALPHA/ (0 . 5*RHOO*VEL.*2*SREF) 

CMALPDT =MALPHDT / (0 . S*RHOO*VEL**2*SREF*CBAR) . (2.VEL/CBAR) 

CZALPDT =ZALPHDT / (0 . 5.RHOO.VEL*.2*SREF) . (2.VEL/CBAR) 

CMQ=MQ/ (0 . 5.RHOO.VEL..2.SREF.CBAR) * (2.VEL/CBAR) 

CZQ=ZQ/ (0 . 5.RHOO.VEL..2.SREF) . (2.VEL/CBAR) 

WRITE (*,200)  CMALPHA, CMALPDT, CMQ 
WRITE (*,201)  CZALPHA, CZALPDT, CZQ 

WRITE(*, ’ (A) ’)  ’  The  control  surface  derivatives  are:’ 

WRITE (4, 200)  CMALPHA, CMALPDT, CMQ 
WRITE (4, 201)  CZALPHA, CZALPDT, CZQ 

WRITE(4, ’ (A) ')  ’  The  control  surface  derivatives  are:’ 

26  CONTINUE 

IF  (NC.EQ.O)  GOTO  7 
DO  7  1=1, NC 

CMDELTA (I) =MDELTA (I) / (0 . 5*RHOO*VEL**2.SREF*CBAR) 

CZDELTA (I) =ZDELTA (I) / (0 . 5*RHOO.VEL*.2*SREF) 

WRITE (* , 202)  I , CMDELTA (I) , I , CZDELTA (I) 

WRITE (4 , 202)  I , CMDELTA (I) , I , CZDELTA (I) 

7  CONTINUE 

IF(NSTAB.EQ.O)  GOTO  8 

WRITE(*, * (/A) ’)  ’  Do  you  want  to  change  any  of  the  derivatives?’ 
READ (*,204)  QCHANGE 
IF (QCHANGE . NE . ’ Y ’ )  GOTO  8 
IF  (NRB.EQ.O)  GOTO  27 
WRITE (*, ’ (A) ’)  ’  CMalpha:’ 

READ(*,*)  CMALPHA 
MALPHA=CMALPH A  *  QBAR  *  SREF  *  CBAR 
WRITE(*, ’ (A) ’)  ’  CZalpha: ’ 

READ(*,*)  CZALPHA 
Z ALPH A=C  Z ALPHA  *  QBAR  *  SREF 
WRITE(*, ’ (A) ’)  ’  CMq: ’ 

READ(*,*)  CMQ 

MQ=CMQ*QBAR*SREF*CBAR**2/ (2*VEL) 

WRITE(*, ’ (A) ’)  ’  CZq: ’ 


READ(*,*)  CZQ 
ZQ=CZQ*SREF*CBAR/ (2*VEL) 

WRITER,* (A)*)  ’  CMalphadt:  * 

READ(*,*)  CMALPDT 

MALPHDT=CMALPDT*Q8AR*SREF*CBAR**2/ (2*VEL) 
WRITE (*, ’ (A) ’)  »  CZalphadt: * 

READ(*,*)  CZALPDT 

ZALPHDT=CZALPDT  *QBAR*SREF*CBAR/ (2*VEL) 


27 


CONTINUE 

IF  (NC.EQ.O)  GOTO  9 
DO  9  1*1, NC 
WRITE (*,205)  I 
READ(*,*)  CMDELTA (I) 

MDELTA (I) =CMDELTA (I) *QBAR*SREF*CBAR 


WRITE (*,206)  I 
READ(*,*)  CZD&TA(I) 

ZDELTA (I) =CZDELTA (I) *QBAR*SREF*CBAR 


9  CONTINUE 

IF  (NRB.EQ.O)  GOTO  28 

WRITE (*, ’ (/A) ’)  *  The  new  der  ivatives  are: 

WRITE (*,200)  CMALPHA , CMALPDT , CMQ 
WRITE (*,201)  CZALPHA, CZALPDT, CZQ 

WRITE(*, ’ (A) ’)  *  The  new  control  surface  derivatives  are: 
WRITE (4, ’ (/A) ’)  ’  The  new  derivatives  are:’ 

WRITE (4, 200)  CMALPHA, CMALPDT, CMQ 
WRITE (4, 201)  CZALPHA .CZALPDT, CZQ 

WRITE (4, ’ (A) ’)  ’  The  new  control  surface  derivatives  are: 
28  CONTINUE 


IF  (NC.EQ.O)  GOTO  10 


DO  10  1=1, NC 

CMDELTA (I) =MDELTA (I) / (0 . 5*RH0O*VEL**2*SREF*CBAR) 
CZDELTA (I) =ZDELTA (I) / (0 . 5*RHOO*VEL*  *2*SREF) 
WRITE (*,202)  I , CMDELTA (I) , I , CZDELTA (I) 

WRITE (4 , 202)  I , CMDELTA (I) , I , CZDELTA (I) 


) 


» 


10  CONTINUE 
8  CONTINUE 

200  FORMAT ( '  CMalphas’.Ell.A, ’ 

201  FORMAT ( ’  CZalpha^Ell.A,  ’ 

202  FORMAT ( ’  CMdelta(Ml,  ’)  =  ’, 

203  FORMAT (//,’  For  a  velocity 
0F8.5, ’slugs/ft**3') 

204  FORMAT (A) 

205  FORMAT ( ’  CMdelta(’,Il, ’) : ’) 

206  FORMAT ( *  CZdel ta(’ ,11, ’) : ’) 


CMa I phadot= ’ , El 1 . 4 , ’  CMq=’,Ell.4) 
CZalphadot=’,Ell.4, ’  CZq=’,Ell.4) 

Ell. 4,’  CZdelta(’,Il, ’)=’ ,E11 .4) 
of ’,F10.4, ’  ft/sec,  and  a  density  of’, 


RETURN 

END 
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